src/gromacs/gmxana/gmx_cluster.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <string.h>
  40 #include <ctype.h>
  41 #include "macros.h"
  42 #include "smalloc.h"
  43 #include "typedefs.h"
  44 #include "statutil.h"
  45 #include "tpxio.h"
  46 #include "string2.h"
  47 #include "vec.h"
  48 #include "macros.h"
  49 #include "index.h"
  50 #include "gmx_random.h"
  51 #include "pbc.h"
  52 #include "rmpbc.h"
  53 #include "xvgr.h"
  54 #include "futil.h"
  55 #include "matio.h"
  56 #include "cmat.h"
  57 #include "do_fit.h"
  58 #include "trnio.h"
  59 #include "viewit.h"
  60 #include "gmx_ana.h"
  61
  62 #include "gromacs/linearalgebra/eigensolver.h"
  63
  64 /* print to two file pointers at once (i.e. stderr and log) */
  65 static inline
  66 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
  67 {
  68     fprintf(fp1, "%s", buf);
  69     fprintf(fp2, "%s", buf);
  70 }
  71
  72 /* just print a prepared buffer to fp1 and fp2 */
  73 static inline
  74 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
  75 {
  76     lo_ffprintf(fp1, fp2, buf);
  77 }
  78
  79 /* prepare buffer with one argument, then print to fp1 and fp2 */
  80 static inline
  81 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
  82 {
  83     sprintf(buf, fmt, arg);
  84     lo_ffprintf(fp1, fp2, buf);
  85 }
  86
  87 /* prepare buffer with one argument, then print to fp1 and fp2 */
  88 static inline
  89 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
  90 {
  91     sprintf(buf, fmt, arg);
  92     lo_ffprintf(fp1, fp2, buf);
  93 }
  94
  95 /* prepare buffer with one argument, then print to fp1 and fp2 */
  96 static inline
  97 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
  98 {
  99     sprintf(buf, fmt, arg);
 100     lo_ffprintf(fp1, fp2, buf);
 101 }
 102
 103 /* prepare buffer with two arguments, then print to fp1 and fp2 */
 104 static inline
 105 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
 106 {
 107     sprintf(buf, fmt, arg1, arg2);
 108     lo_ffprintf(fp1, fp2, buf);
 109 }
 110
 111 /* prepare buffer with two arguments, then print to fp1 and fp2 */
 112 static inline
 113 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
 114 {
 115     sprintf(buf, fmt, arg1, arg2);
 116     lo_ffprintf(fp1, fp2, buf);
 117 }
 118
 119 /* prepare buffer with two arguments, then print to fp1 and fp2 */
 120 static inline
 121 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
 122 {
 123     sprintf(buf, fmt, arg1, arg2);
 124     lo_ffprintf(fp1, fp2, buf);
 125 }
 126
 127 typedef struct {
 128     int  ncl;
 129     int *cl;
 130 } t_clusters;
 131
 132 typedef struct {
 133     int  nr;
 134     int *nb;
 135 } t_nnb;
 136
 137 void pr_energy(FILE *fp, real e)
 138 {
 139     fprintf(fp, "Energy: %8.4f\n", e);
 140 }
 141
 142 void cp_index(int nn, int from[], int to[])
 143 {
 144     int i;
 145
 146     for (i = 0; (i < nn); i++)
 147     {
 148         to[i] = from[i];
 149     }
 150 }
 151
 152 void mc_optimize(FILE *log, t_mat *m, real *time,
 153                  int maxiter, int nrandom,
 154                  int seed, real kT,
 155                  const char *conv, output_env_t oenv)
 156 {
 157     FILE      *fp = NULL;
 158     real       ecur, enext, emin, prob, enorm;
 159     int        i, j, iswap, jswap, nn, nuphill = 0;
 160     gmx_rng_t  rng;
 161     t_mat     *minimum;
 162
 163     if (m->n1 != m->nn)
 164     {
 165         fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
 166         return;
 167     }
 168     printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
 169     printf("by reordering the frames to minimize the path between the two structures\n");
 170     printf("that have the largest pairwise RMSD.\n");
 171
 172     iswap = jswap = -1;
 173     enorm = m->mat[0][0];
 174     for (i = 0; (i < m->n1); i++)
 175     {
 176         for (j = 0; (j < m->nn); j++)
 177         {
 178             if (m->mat[i][j] > enorm)
 179             {
 180                 enorm   = m->mat[i][j];
 181                 iswap   = i;
 182                 jswap   = j;
 183             }
 184         }
 185     }
 186     if ((iswap == -1) || (jswap == -1))
 187     {
 188         fprintf(stderr, "Matrix contains identical values in all fields\n");
 189         return;
 190     }
 191     swap_rows(m, 0, iswap);
 192     swap_rows(m, m->n1-1, jswap);
 193     emin = ecur = mat_energy(m);
 194     printf("Largest distance %g between %d and %d. Energy: %g.\n",
 195            enorm, iswap, jswap, emin);
 196
 197     rng = gmx_rng_init(seed);
 198     nn  = m->nn;
 199
 200     /* Initiate and store global minimum */
 201     minimum     = init_mat(nn, m->b1D);
 202     minimum->nn = nn;
 203     copy_t_mat(minimum, m);
 204
 205     if (NULL != conv)
 206     {
 207         fp = xvgropen(conv, "Convergence of the MC optimization",
 208                       "Energy", "Step", oenv);
 209     }
 210     for (i = 0; (i < maxiter); i++)
 211     {
 212         /* Generate new swapping candidates */
 213         do
 214         {
 215             iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
 216             jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
 217         }
 218         while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
 219
 220         /* Apply swap and compute energy */
 221         swap_rows(m, iswap, jswap);
 222         enext = mat_energy(m);
 223
 224         /* Compute probability */
 225         prob = 0;
 226         if ((enext < ecur) || (i < nrandom))
 227         {
 228             prob = 1;
 229             if (enext < emin)
 230             {
 231                 /* Store global minimum */
 232                 copy_t_mat(minimum, m);
 233                 emin = enext;
 234             }
 235         }
 236         else if (kT > 0)
 237         {
 238             /* Try Monte Carlo step */
 239             prob = exp(-(enext-ecur)/(enorm*kT));
 240         }
 241
 242         if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
 243         {
 244             if (enext > ecur)
 245             {
 246                 nuphill++;
 247             }
 248
 249             fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
 250                     i, iswap, jswap, enext, prob);
 251             if (NULL != fp)
 252             {
 253                 fprintf(fp, "%6d  %10g\n", i, enext);
 254             }
 255             ecur = enext;
 256         }
 257         else
 258         {
 259             swap_rows(m, jswap, iswap);
 260         }
 261     }
 262     fprintf(log, "%d uphill steps were taken during optimization\n",
 263             nuphill);
 264
 265     /* Now swap the matrix to get it into global minimum mode */
 266     copy_t_mat(m, minimum);
 267
 268     fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
 269     fprintf(log, "Global minimum energy %g\n", mat_energy(m));
 270     fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
 271     for (i = 0; (i < m->nn); i++)
 272     {
 273         fprintf(log, "%10g  %5d  %10g\n",
 274                 time[m->m_ind[i]],
 275                 m->m_ind[i],
 276                 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
 277     }
 278
 279     if (NULL != fp)
 280     {
 281         fclose(fp);
 282     }
 283 }
 284
 285 static void calc_dist(int nind, rvec x[], real **d)
 286 {
 287     int      i, j;
 288     real    *xi;
 289     rvec     dx;
 290
 291     for (i = 0; (i < nind-1); i++)
 292     {
 293         xi = x[i];
 294         for (j = i+1; (j < nind); j++)
 295         {
 296             /* Should use pbc_dx when analysing multiple molecueles,
 297              * but the box is not stored for every frame.
 298              */
 299             rvec_sub(xi, x[j], dx);
 300             d[i][j] = norm(dx);
 301         }
 302     }
 303 }
 304
 305 static real rms_dist(int isize, real **d, real **d_r)
 306 {
 307     int  i, j;
 308     real r, r2;
 309
 310     r2 = 0.0;
 311     for (i = 0; (i < isize-1); i++)
 312     {
 313         for (j = i+1; (j < isize); j++)
 314         {
 315             r   = d[i][j]-d_r[i][j];
 316             r2 += r*r;
 317         }
 318     }
 319     r2 /= (isize*(isize-1))/2;
 320
 321     return sqrt(r2);
 322 }
 323
 324 static int rms_dist_comp(const void *a, const void *b)
 325 {
 326     t_dist *da, *db;
 327
 328     da = (t_dist *)a;
 329     db = (t_dist *)b;
 330
 331     if (da->dist - db->dist < 0)
 332     {
 333         return -1;
 334     }
 335     else if (da->dist - db->dist > 0)
 336     {
 337         return 1;
 338     }
 339     return 0;
 340 }
 341
 342 static int clust_id_comp(const void *a, const void *b)
 343 {
 344     t_clustid *da, *db;
 345
 346     da = (t_clustid *)a;
 347     db = (t_clustid *)b;
 348
 349     return da->clust - db->clust;
 350 }
 351
 352 static int nrnb_comp(const void *a, const void *b)
 353 {
 354     t_nnb *da, *db;
 355
 356     da = (t_nnb *)a;
 357     db = (t_nnb *)b;
 358
 359     /* return the b-a, we want highest first */
 360     return db->nr - da->nr;
 361 }
 362
 363 void gather(t_mat *m, real cutoff, t_clusters *clust)
 364 {
 365     t_clustid *c;
 366     t_dist    *d;
 367     int        i, j, k, nn, cid, n1, diff;
 368     gmx_bool   bChange;
 369
 370     /* First we sort the entries in the RMSD matrix */
 371     n1 = m->nn;
 372     nn = ((n1-1)*n1)/2;
 373     snew(d, nn);
 374     for (i = k = 0; (i < n1); i++)
 375     {
 376         for (j = i+1; (j < n1); j++, k++)
 377         {
 378             d[k].i    = i;
 379             d[k].j    = j;
 380             d[k].dist = m->mat[i][j];
 381         }
 382     }
 383     if (k != nn)
 384     {
 385         gmx_incons("gather algortihm");
 386     }
 387     qsort(d, nn, sizeof(d[0]), rms_dist_comp);
 388
 389     /* Now we make a cluster index for all of the conformations */
 390     c = new_clustid(n1);
 391
 392     /* Now we check the closest structures, and equalize their cluster numbers */
 393     fprintf(stderr, "Linking structures ");
 394     do
 395     {
 396         fprintf(stderr, "*");
 397         bChange = FALSE;
 398         for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
 399         {
 400             diff = c[d[k].j].clust - c[d[k].i].clust;
 401             if (diff)
 402             {
 403                 bChange = TRUE;
 404                 if (diff > 0)
 405                 {
 406                     c[d[k].j].clust = c[d[k].i].clust;
 407                 }
 408                 else
 409                 {
 410                     c[d[k].i].clust = c[d[k].j].clust;
 411                 }
 412             }
 413         }
 414     }
 415     while (bChange);
 416     fprintf(stderr, "\nSorting and renumbering clusters\n");
 417     /* Sort on cluster number */
 418     qsort(c, n1, sizeof(c[0]), clust_id_comp);
 419
 420     /* Renumber clusters */
 421     cid = 1;
 422     for (k = 1; k < n1; k++)
 423     {
 424         if (c[k].clust != c[k-1].clust)
 425         {
 426             c[k-1].clust = cid;
 427             cid++;
 428         }
 429         else
 430         {
 431             c[k-1].clust = cid;
 432         }
 433     }
 434     c[k-1].clust = cid;
 435     if (debug)
 436     {
 437         for (k = 0; (k < n1); k++)
 438         {
 439             fprintf(debug, "Cluster index for conformation %d: %d\n",
 440                     c[k].conf, c[k].clust);
 441         }
 442     }
 443     clust->ncl = cid;
 444     for (k = 0; k < n1; k++)
 445     {
 446         clust->cl[c[k].conf] = c[k].clust;
 447     }
 448
 449     sfree(c);
 450     sfree(d);
 451 }
 452
 453 gmx_bool jp_same(int **nnb, int i, int j, int P)
 454 {
 455     gmx_bool bIn;
 456     int      k, ii, jj, pp;
 457
 458     bIn = FALSE;
 459     for (k = 0; nnb[i][k] >= 0; k++)
 460     {
 461         bIn = bIn || (nnb[i][k] == j);
 462     }
 463     if (!bIn)
 464     {
 465         return FALSE;
 466     }
 467
 468     bIn = FALSE;
 469     for (k = 0; nnb[j][k] >= 0; k++)
 470     {
 471         bIn = bIn || (nnb[j][k] == i);
 472     }
 473     if (!bIn)
 474     {
 475         return FALSE;
 476     }
 477
 478     pp = 0;
 479     for (ii = 0; nnb[i][ii] >= 0; ii++)
 480     {
 481         for (jj = 0; nnb[j][jj] >= 0; jj++)
 482         {
 483             if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
 484             {
 485                 pp++;
 486             }
 487         }
 488     }
 489
 490     return (pp >= P);
 491 }
 492
 493 static void jarvis_patrick(int n1, real **mat, int M, int P,
 494                            real rmsdcut, t_clusters *clust)
 495 {
 496     t_dist     *row;
 497     t_clustid  *c;
 498     int       **nnb;
 499     int         i, j, k, cid, diff, max;
 500     gmx_bool    bChange;
 501     real      **mcpy = NULL;
 502
 503     if (rmsdcut < 0)
 504     {
 505         rmsdcut = 10000;
 506     }
 507
 508     /* First we sort the entries in the RMSD matrix row by row.
 509      * This gives us the nearest neighbor list.
 510      */
 511     snew(nnb, n1);
 512     snew(row, n1);
 513     for (i = 0; (i < n1); i++)
 514     {
 515         for (j = 0; (j < n1); j++)
 516         {
 517             row[j].j    = j;
 518             row[j].dist = mat[i][j];
 519         }
 520         qsort(row, n1, sizeof(row[0]), rms_dist_comp);
 521         if (M > 0)
 522         {
 523             /* Put the M nearest neighbors in the list */
 524             snew(nnb[i], M+1);
 525             for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
 526             {
 527                 if (row[j].j  != i)
 528                 {
 529                     nnb[i][k]  = row[j].j;
 530                     k++;
 531                 }
 532             }
 533             nnb[i][k] = -1;
 534         }
 535         else
 536         {
 537             /* Put all neighbors nearer than rmsdcut in the list */
 538             max = 0;
 539             k   = 0;
 540             for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
 541             {
 542                 if (row[j].j != i)
 543                 {
 544                     if (k >= max)
 545                     {
 546                         max += 10;
 547                         srenew(nnb[i], max);
 548                     }
 549                     nnb[i][k] = row[j].j;
 550                     k++;
 551                 }
 552             }
 553             if (k == max)
 554             {
 555                 srenew(nnb[i], max+1);
 556             }
 557             nnb[i][k] = -1;
 558         }
 559     }
 560     sfree(row);
 561     if (debug)
 562     {
 563         fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
 564         for (i = 0; (i < n1); i++)
 565         {
 566             fprintf(debug, "i:%5d nbs:", i);
 567             for (j = 0; nnb[i][j] >= 0; j++)
 568             {
 569                 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
 570             }
 571             fprintf(debug, "\n");
 572         }
 573     }
 574
 575     c = new_clustid(n1);
 576     fprintf(stderr, "Linking structures ");
 577     /* Use mcpy for temporary storage of booleans */
 578     mcpy = mk_matrix(n1, n1, FALSE);
 579     for (i = 0; i < n1; i++)
 580     {
 581         for (j = i+1; j < n1; j++)
 582         {
 583             mcpy[i][j] = jp_same(nnb, i, j, P);
 584         }
 585     }
 586     do
 587     {
 588         fprintf(stderr, "*");
 589         bChange = FALSE;
 590         for (i = 0; i < n1; i++)
 591         {
 592             for (j = i+1; j < n1; j++)
 593             {
 594                 if (mcpy[i][j])
 595                 {
 596                     diff = c[j].clust - c[i].clust;
 597                     if (diff)
 598                     {
 599                         bChange = TRUE;
 600                         if (diff > 0)
 601                         {
 602                             c[j].clust = c[i].clust;
 603                         }
 604                         else
 605                         {
 606                             c[i].clust = c[j].clust;
 607                         }
 608                     }
 609                 }
 610             }
 611         }
 612     }
 613     while (bChange);
 614
 615     fprintf(stderr, "\nSorting and renumbering clusters\n");
 616     /* Sort on cluster number */
 617     qsort(c, n1, sizeof(c[0]), clust_id_comp);
 618
 619     /* Renumber clusters */
 620     cid = 1;
 621     for (k = 1; k < n1; k++)
 622     {
 623         if (c[k].clust != c[k-1].clust)
 624         {
 625             c[k-1].clust = cid;
 626             cid++;
 627         }
 628         else
 629         {
 630             c[k-1].clust = cid;
 631         }
 632     }
 633     c[k-1].clust = cid;
 634     clust->ncl   = cid;
 635     for (k = 0; k < n1; k++)
 636     {
 637         clust->cl[c[k].conf] = c[k].clust;
 638     }
 639     if (debug)
 640     {
 641         for (k = 0; (k < n1); k++)
 642         {
 643             fprintf(debug, "Cluster index for conformation %d: %d\n",
 644                     c[k].conf, c[k].clust);
 645         }
 646     }
 647
 648 /* Again, I don't see the point in this... (AF) */
 649 /*   for(i=0; (i<n1); i++) { */
 650 /*     for(j=0; (j<n1); j++) */
 651 /*       mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
 652 /*   } */
 653 /*   for(i=0; (i<n1); i++) { */
 654 /*     for(j=0; (j<n1); j++) */
 655 /*       mat[i][j] = mcpy[i][j]; */
 656 /*   } */
 657     done_matrix(n1, &mcpy);
 658
 659     sfree(c);
 660     for (i = 0; (i < n1); i++)
 661     {
 662         sfree(nnb[i]);
 663     }
 664     sfree(nnb);
 665 }
 666
 667 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
 668 {
 669     int i, j;
 670
 671     /* dump neighbor list */
 672     fprintf(fp, "%s", title);
 673     for (i = 0; (i < n1); i++)
 674     {
 675         fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
 676         for (j = 0; j < nnb[i].nr; j++)
 677         {
 678             fprintf(fp, "%5d", nnb[i].nb[j]);
 679         }
 680         fprintf(fp, "\n");
 681     }
 682 }
 683
 684 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
 685 {
 686     t_dist *row;
 687     t_nnb  *nnb;
 688     int     i, j, k, j1, max;
 689
 690     /* Put all neighbors nearer than rmsdcut in the list */
 691     fprintf(stderr, "Making list of neighbors within cutoff ");
 692     snew(nnb, n1);
 693     snew(row, n1);
 694     for (i = 0; (i < n1); i++)
 695     {
 696         max = 0;
 697         k   = 0;
 698         /* put all neighbors within cut-off in list */
 699         for (j = 0; j < n1; j++)
 700         {
 701             if (mat[i][j] < rmsdcut)
 702             {
 703                 if (k >= max)
 704                 {
 705                     max += 10;
 706                     srenew(nnb[i].nb, max);
 707                 }
 708                 nnb[i].nb[k] = j;
 709                 k++;
 710             }
 711         }
 712         /* store nr of neighbors, we'll need that */
 713         nnb[i].nr = k;
 714         if (i%(1+n1/100) == 0)
 715         {
 716             fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
 717         }
 718     }
 719     fprintf(stderr, "%3d%%\n", 100);
 720     sfree(row);
 721
 722     /* sort neighbor list on number of neighbors, largest first */
 723     qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
 724
 725     if (debug)
 726     {
 727         dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
 728     }
 729
 730     /* turn first structure with all its neighbors (largest) into cluster
 731        remove them from pool of structures and repeat for all remaining */
 732     fprintf(stderr, "Finding clusters %4d", 0);
 733     /* cluster id's start at 1: */
 734     k = 1;
 735     while (nnb[0].nr)
 736     {
 737         /* set cluster id (k) for first item in neighborlist */
 738         for (j = 0; j < nnb[0].nr; j++)
 739         {
 740             clust->cl[nnb[0].nb[j]] = k;
 741         }
 742         /* mark as done */
 743         nnb[0].nr = 0;
 744         sfree(nnb[0].nb);
 745
 746         /* adjust number of neighbors for others, taking removals into account: */
 747         for (i = 1; i < n1 && nnb[i].nr; i++)
 748         {
 749             j1 = 0;
 750             for (j = 0; j < nnb[i].nr; j++)
 751             {
 752                 /* if this neighbor wasn't removed */
 753                 if (clust->cl[nnb[i].nb[j]] == 0)
 754                 {
 755                     /* shift the rest (j1<=j) */
 756                     nnb[i].nb[j1] = nnb[i].nb[j];
 757                     /* next */
 758                     j1++;
 759                 }
 760             }
 761             /* now j1 is the new number of neighbors */
 762             nnb[i].nr = j1;
 763         }
 764         /* sort again on nnb[].nr, because we have new # neighbors: */
 765         /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
 766         qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
 767
 768         fprintf(stderr, "\b\b\b\b%4d", k);
 769         /* new cluster id */
 770         k++;
 771     }
 772     fprintf(stderr, "\n");
 773     sfree(nnb);
 774     if (debug)
 775     {
 776         fprintf(debug, "Clusters (%d):\n", k);
 777         for (i = 0; i < n1; i++)
 778         {
 779             fprintf(debug, " %3d", clust->cl[i]);
 780         }
 781         fprintf(debug, "\n");
 782     }
 783
 784     clust->ncl = k-1;
 785 }
 786
 787 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
 788                       int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
 789 {
 790     rvec       **xx, *x;
 791     matrix       box;
 792     real         t;
 793     int          i, i0, j, max_nf;
 794     int          natom;
 795     t_trxstatus *status;
 796
 797
 798     max_nf = 0;
 799     xx     = NULL;
 800     *time  = NULL;
 801     natom  = read_first_x(oenv, &status, fn, &t, &x, box);
 802     i      = 0;
 803     i0     = 0;
 804     do
 805     {
 806         if (bPBC)
 807         {
 808             gmx_rmpbc(gpbc, natom, box, x);
 809         }
 810         if (i0 >= max_nf)
 811         {
 812             max_nf += 10;
 813             srenew(xx, max_nf);
 814             srenew(*time, max_nf);
 815         }
 816         if ((i % skip) == 0)
 817         {
 818             snew(xx[i0], isize);
 819             /* Store only the interesting atoms */
 820             for (j = 0; (j < isize); j++)
 821             {
 822                 copy_rvec(x[index[j]], xx[i0][j]);
 823             }
 824             (*time)[i0] = t;
 825             i0++;
 826         }
 827         i++;
 828     }
 829     while (read_next_x(oenv, status, &t, x, box));
 830     fprintf(stderr, "Allocated %lu bytes for frames\n",
 831             (unsigned long) (max_nf*isize*sizeof(**xx)));
 832     fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
 833     *nframe = i0;
 834     sfree(x);
 835
 836     return xx;
 837 }
 838
 839 static int plot_clusters(int nf, real **mat, t_clusters *clust,
 840                          int minstruct)
 841 {
 842     int  i, j, ncluster, ci;
 843     int *cl_id, *nstruct, *strind;
 844
 845     snew(cl_id, nf);
 846     snew(nstruct, nf);
 847     snew(strind, nf);
 848     for (i = 0; i < nf; i++)
 849     {
 850         strind[i] = 0;
 851         cl_id[i]  = clust->cl[i];
 852         nstruct[cl_id[i]]++;
 853     }
 854     ncluster = 0;
 855     for (i = 0; i < nf; i++)
 856     {
 857         if (nstruct[i] >= minstruct)
 858         {
 859             ncluster++;
 860             for (j = 0; (j < nf); j++)
 861             {
 862                 if (cl_id[j] == i)
 863                 {
 864                     strind[j] = ncluster;
 865                 }
 866             }
 867         }
 868     }
 869     ncluster++;
 870     fprintf(stderr, "There are %d clusters with at least %d conformations\n",
 871             ncluster, minstruct);
 872
 873     for (i = 0; (i < nf); i++)
 874     {
 875         ci = cl_id[i];
 876         for (j = 0; j < i; j++)
 877         {
 878             if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
 879             {
 880                 /* color different clusters with different colors, as long as
 881                    we don't run out of colors */
 882                 mat[i][j] = strind[i];
 883             }
 884             else
 885             {
 886                 mat[i][j] = 0;
 887             }
 888         }
 889     }
 890     sfree(strind);
 891     sfree(nstruct);
 892     sfree(cl_id);
 893
 894     return ncluster;
 895 }
 896
 897 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
 898 {
 899     int i, j, v;
 900
 901     for (i = 0; i < nf; i++)
 902     {
 903         for (j = 0; j < i; j++)
 904         {
 905             if (clust->cl[i] == clust->cl[j])
 906             {
 907                 mat[i][j] = val;
 908             }
 909             else
 910             {
 911                 mat[i][j] = 0;
 912             }
 913         }
 914     }
 915 }
 916
 917 static char *parse_filename(const char *fn, int maxnr)
 918 {
 919     int         i;
 920     char       *fnout;
 921     const char *ext;
 922     char        buf[STRLEN];
 923
 924     if (strchr(fn, '%'))
 925     {
 926         gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
 927     }
 928     /* number of digits needed in numbering */
 929     i = (int)(log(maxnr)/log(10)) + 1;
 930     /* split fn and ext */
 931     ext = strrchr(fn, '.');
 932     if (!ext)
 933     {
 934         gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
 935     }
 936     ext++;
 937     /* insert e.g. '%03d' between fn and ext */
 938     sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
 939     snew(fnout, strlen(buf)+1);
 940     strcpy(fnout, buf);
 941
 942     return fnout;
 943 }
 944
 945 static void ana_trans(t_clusters *clust, int nf,
 946                       const char *transfn, const char *ntransfn, FILE *log,
 947                       t_rgb rlo, t_rgb rhi, const output_env_t oenv)
 948 {
 949     FILE  *fp;
 950     real **trans, *axis;
 951     int   *ntrans;
 952     int    i, ntranst, maxtrans;
 953     char   buf[STRLEN];
 954
 955     snew(ntrans, clust->ncl);
 956     snew(trans, clust->ncl);
 957     snew(axis, clust->ncl);
 958     for (i = 0; i < clust->ncl; i++)
 959     {
 960         axis[i] = i+1;
 961         snew(trans[i], clust->ncl);
 962     }
 963     ntranst  = 0;
 964     maxtrans = 0;
 965     for (i = 1; i < nf; i++)
 966     {
 967         if (clust->cl[i] != clust->cl[i-1])
 968         {
 969             ntranst++;
 970             ntrans[clust->cl[i-1]-1]++;
 971             ntrans[clust->cl[i]-1]++;
 972             trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
 973             maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
 974         }
 975     }
 976     ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
 977                 "max %d between two specific clusters\n", ntranst, maxtrans);
 978     if (transfn)
 979     {
 980         fp = ffopen(transfn, "w");
 981         i  = min(maxtrans+1, 80);
 982         write_xpm(fp, 0, "Cluster Transitions", "# transitions",
 983                   "from cluster", "to cluster",
 984                   clust->ncl, clust->ncl, axis, axis, trans,
 985                   0, maxtrans, rlo, rhi, &i);
 986         ffclose(fp);
 987     }
 988     if (ntransfn)
 989     {
 990         fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
 991                       oenv);
 992         for (i = 0; i < clust->ncl; i++)
 993         {
 994             fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
 995         }
 996         ffclose(fp);
 997     }
 998     sfree(ntrans);
 999     for (i = 0; i < clust->ncl; i++)
1000     {
1001         sfree(trans[i]);
1002     }
1003     sfree(trans);
1004     sfree(axis);
1005 }
1006
1007 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1008                              int natom, t_atoms *atoms, rvec *xtps,
1009                              real *mass, rvec **xx, real *time,
1010                              int ifsize, atom_id *fitidx,
1011                              int iosize, atom_id *outidx,
1012                              const char *trxfn, const char *sizefn,
1013                              const char *transfn, const char *ntransfn,
1014                              const char *clustidfn, gmx_bool bAverage,
1015                              int write_ncl, int write_nst, real rmsmin,
1016                              gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1017                              const output_env_t oenv)
1018 {
1019     FILE        *fp = NULL;
1020     char         buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1021     t_trxstatus *trxout  = NULL;
1022     t_trxstatus *trxsout = NULL;
1023     int          i, i1, cl, nstr, *structure, first = 0, midstr;
1024     gmx_bool    *bWrite = NULL;
1025     real         r, clrmsd, midrmsd;
1026     rvec        *xav = NULL;
1027     matrix       zerobox;
1028
1029     clear_mat(zerobox);
1030
1031     ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1032     trxsfn = NULL;
1033     if (trxfn)
1034     {
1035         /* do we write all structures? */
1036         if (write_ncl)
1037         {
1038             trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1039             snew(bWrite, nf);
1040         }
1041         ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1042                     bAverage ? "average" : "middle", trxfn);
1043         if (write_ncl)
1044         {
1045             /* find out what we want to tell the user:
1046                Writing [all structures|structures with rmsd > %g] for
1047                {all|first %d} clusters {with more than %d structures} to %s     */
1048             if (rmsmin > 0.0)
1049             {
1050                 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1051             }
1052             else
1053             {
1054                 sprintf(buf1, "all structures");
1055             }
1056             buf2[0] = buf3[0] = '\0';
1057             if (write_ncl >= clust->ncl)
1058             {
1059                 if (write_nst == 0)
1060                 {
1061                     sprintf(buf2, "all ");
1062                 }
1063             }
1064             else
1065             {
1066                 sprintf(buf2, "the first %d ", write_ncl);
1067             }
1068             if (write_nst)
1069             {
1070                 sprintf(buf3, " with more than %d structures", write_nst);
1071             }
1072             sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1073             ffprintf(stderr, log, buf);
1074         }
1075
1076         /* Prepare a reference structure for the orientation of the clusters  */
1077         if (bFit)
1078         {
1079             reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1080         }
1081         trxout = open_trx(trxfn, "w");
1082         /* Calculate the average structure in each cluster,               *
1083          * all structures are fitted to the first struture of the cluster */
1084         snew(xav, natom);
1085     }
1086
1087     if (transfn || ntransfn)
1088     {
1089         ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1090     }
1091
1092     if (clustidfn)
1093     {
1094         fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1095         fprintf(fp, "@    s0 symbol 2\n");
1096         fprintf(fp, "@    s0 symbol size 0.2\n");
1097         fprintf(fp, "@    s0 linestyle 0\n");
1098         for (i = 0; i < nf; i++)
1099         {
1100             fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1101         }
1102         ffclose(fp);
1103     }
1104     if (sizefn)
1105     {
1106         fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1107         fprintf(fp, "@g%d type %s\n", 0, "bar");
1108     }
1109     snew(structure, nf);
1110     fprintf(log, "\n%3s | %3s  %4s | %6s %4s | cluster members\n",
1111             "cl.", "#st", "rmsd", "middle", "rmsd");
1112     for (cl = 1; cl <= clust->ncl; cl++)
1113     {
1114         /* prepare structures (fit, middle, average) */
1115         if (xav)
1116         {
1117             for (i = 0; i < natom; i++)
1118             {
1119                 clear_rvec(xav[i]);
1120             }
1121         }
1122         nstr = 0;
1123         for (i1 = 0; i1 < nf; i1++)
1124         {
1125             if (clust->cl[i1] == cl)
1126             {
1127                 structure[nstr] = i1;
1128                 nstr++;
1129                 if (trxfn && (bAverage || write_ncl) )
1130                 {
1131                     if (bFit)
1132                     {
1133                         reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1134                     }
1135                     if (nstr == 1)
1136                     {
1137                         first = i1;
1138                     }
1139                     else if (bFit)
1140                     {
1141                         do_fit(natom, mass, xx[first], xx[i1]);
1142                     }
1143                     if (xav)
1144                     {
1145                         for (i = 0; i < natom; i++)
1146                         {
1147                             rvec_inc(xav[i], xx[i1][i]);
1148                         }
1149                     }
1150                 }
1151             }
1152         }
1153         if (sizefn)
1154         {
1155             fprintf(fp, "%8d %8d\n", cl, nstr);
1156         }
1157         clrmsd  = 0;
1158         midstr  = 0;
1159         midrmsd = 10000;
1160         for (i1 = 0; i1 < nstr; i1++)
1161         {
1162             r = 0;
1163             if (nstr > 1)
1164             {
1165                 for (i = 0; i < nstr; i++)
1166                 {
1167                     if (i < i1)
1168                     {
1169                         r += rmsd[structure[i]][structure[i1]];
1170                     }
1171                     else
1172                     {
1173                         r += rmsd[structure[i1]][structure[i]];
1174                     }
1175                 }
1176                 r /= (nstr - 1);
1177             }
1178             if (r < midrmsd)
1179             {
1180                 midstr  = structure[i1];
1181                 midrmsd = r;
1182             }
1183             clrmsd += r;
1184         }
1185         clrmsd /= nstr;
1186
1187         /* dump cluster info to logfile */
1188         if (nstr > 1)
1189         {
1190             sprintf(buf1, "%6.3f", clrmsd);
1191             if (buf1[0] == '0')
1192             {
1193                 buf1[0] = ' ';
1194             }
1195             sprintf(buf2, "%5.3f", midrmsd);
1196             if (buf2[0] == '0')
1197             {
1198                 buf2[0] = ' ';
1199             }
1200         }
1201         else
1202         {
1203             sprintf(buf1, "%5s", "");
1204             sprintf(buf2, "%5s", "");
1205         }
1206         fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1207         for (i = 0; i < nstr; i++)
1208         {
1209             if ((i % 7 == 0) && i)
1210             {
1211                 sprintf(buf, "\n%3s | %3s  %4s | %6s %4s |", "", "", "", "", "");
1212             }
1213             else
1214             {
1215                 buf[0] = '\0';
1216             }
1217             i1 = structure[i];
1218             fprintf(log, "%s %6g", buf, time[i1]);
1219         }
1220         fprintf(log, "\n");
1221
1222         /* write structures to trajectory file(s) */
1223         if (trxfn)
1224         {
1225             if (write_ncl)
1226             {
1227                 for (i = 0; i < nstr; i++)
1228                 {
1229                     bWrite[i] = FALSE;
1230                 }
1231             }
1232             if (cl < write_ncl+1 && nstr > write_nst)
1233             {
1234                 /* Dump all structures for this cluster */
1235                 /* generate numbered filename (there is a %d in trxfn!) */
1236                 sprintf(buf, trxsfn, cl);
1237                 trxsout = open_trx(buf, "w");
1238                 for (i = 0; i < nstr; i++)
1239                 {
1240                     bWrite[i] = TRUE;
1241                     if (rmsmin > 0.0)
1242                     {
1243                         for (i1 = 0; i1 < i && bWrite[i]; i1++)
1244                         {
1245                             if (bWrite[i1])
1246                             {
1247                                 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1248                             }
1249                         }
1250                     }
1251                     if (bWrite[i])
1252                     {
1253                         write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1254                                   xx[structure[i]], NULL, NULL);
1255                     }
1256                 }
1257                 close_trx(trxsout);
1258             }
1259             /* Dump the average structure for this cluster */
1260             if (bAverage)
1261             {
1262                 for (i = 0; i < natom; i++)
1263                 {
1264                     svmul(1.0/nstr, xav[i], xav[i]);
1265                 }
1266             }
1267             else
1268             {
1269                 for (i = 0; i < natom; i++)
1270                 {
1271                     copy_rvec(xx[midstr][i], xav[i]);
1272                 }
1273                 if (bFit)
1274                 {
1275                     reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1276                 }
1277             }
1278             if (bFit)
1279             {
1280                 do_fit(natom, mass, xtps, xav);
1281             }
1282             r = cl;
1283             write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1284         }
1285     }
1286     /* clean up */
1287     if (trxfn)
1288     {
1289         close_trx(trxout);
1290         sfree(xav);
1291         if (write_ncl)
1292         {
1293             sfree(bWrite);
1294         }
1295     }
1296     sfree(structure);
1297     if (trxsfn)
1298     {
1299         sfree(trxsfn);
1300     }
1301 }
1302
1303 static void convert_mat(t_matrix *mat, t_mat *rms)
1304 {
1305     int i, j;
1306
1307     rms->n1 = mat->nx;
1308     matrix2real(mat, rms->mat);
1309     /* free input xpm matrix data */
1310     for (i = 0; i < mat->nx; i++)
1311     {
1312         sfree(mat->matrix[i]);
1313     }
1314     sfree(mat->matrix);
1315
1316     for (i = 0; i < mat->nx; i++)
1317     {
1318         for (j = i; j < mat->nx; j++)
1319         {
1320             rms->sumrms += rms->mat[i][j];
1321             rms->maxrms  = max(rms->maxrms, rms->mat[i][j]);
1322             if (j != i)
1323             {
1324                 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1325             }
1326         }
1327     }
1328     rms->nn = mat->nx;
1329 }
1330
1331 int gmx_cluster(int argc, char *argv[])
1332 {
1333     const char        *desc[] = {
1334         "[TT]g_cluster[tt] can cluster structures using several different methods.",
1335         "Distances between structures can be determined from a trajectory",
1336         "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1337         "RMS deviation after fitting or RMS deviation of atom-pair distances",
1338         "can be used to define the distance between structures.[PAR]",
1339
1340         "single linkage: add a structure to a cluster when its distance to any",
1341         "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1342
1343         "Jarvis Patrick: add a structure to a cluster when this structure",
1344         "and a structure in the cluster have each other as neighbors and",
1345         "they have a least [TT]P[tt] neighbors in common. The neighbors",
1346         "of a structure are the M closest structures or all structures within",
1347         "[TT]cutoff[tt].[PAR]",
1348
1349         "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1350         "the order of the frames is using the smallest possible increments.",
1351         "With this it is possible to make a smooth animation going from one",
1352         "structure to another with the largest possible (e.g.) RMSD between",
1353         "them, however the intermediate steps should be as small as possible.",
1354         "Applications could be to visualize a potential of mean force",
1355         "ensemble of simulations or a pulling simulation. Obviously the user",
1356         "has to prepare the trajectory well (e.g. by not superimposing frames).",
1357         "The final result can be inspect visually by looking at the matrix",
1358         "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1359
1360         "diagonalization: diagonalize the RMSD matrix.[PAR]",
1361
1362         "gromos: use algorithm as described in Daura [IT]et al.[it]",
1363         "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1364         "Count number of neighbors using cut-off, take structure with",
1365         "largest number of neighbors with all its neighbors as cluster",
1366         "and eliminate it from the pool of clusters. Repeat for remaining",
1367         "structures in pool.[PAR]",
1368
1369         "When the clustering algorithm assigns each structure to exactly one",
1370         "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1371         "file is supplied, the structure with",
1372         "the smallest average distance to the others or the average structure",
1373         "or all structures for each cluster will be written to a trajectory",
1374         "file. When writing all structures, separate numbered files are made",
1375         "for each cluster.[PAR]",
1376
1377         "Two output files are always written:[BR]",
1378         "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1379         "and a graphical depiction of the clusters in the lower right half",
1380         "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1381         "when two structures are in the same cluster.",
1382         "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1383         "cluster.[BR]",
1384         "[TT]-g[tt] writes information on the options used and a detailed list",
1385         "of all clusters and their members.[PAR]",
1386
1387         "Additionally, a number of optional output files can be written:[BR]",
1388         "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1389         "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1390         "diagonalization.[BR]",
1391         "[TT]-sz[tt] writes the cluster sizes.[BR]",
1392         "[TT]-tr[tt] writes a matrix of the number transitions between",
1393         "cluster pairs.[BR]",
1394         "[TT]-ntr[tt] writes the total number of transitions to or from",
1395         "each cluster.[BR]",
1396         "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1397         "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1398         "structure of each cluster or writes numbered files with cluster members",
1399         "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1400         "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1401         "structure with the smallest average RMSD from all other structures",
1402         "of the cluster.[BR]",
1403     };
1404
1405     FILE              *fp, *log;
1406     int                nf, i, i1, i2, j;
1407     gmx_large_int_t    nrms = 0;
1408
1409     matrix             box;
1410     rvec              *xtps, *usextps, *x1, **xx = NULL;
1411     const char        *fn, *trx_out_fn;
1412     t_clusters         clust;
1413     t_mat             *rms, *orig = NULL;
1414     real              *eigenvalues;
1415     t_topology         top;
1416     int                ePBC;
1417     t_atoms            useatoms;
1418     t_matrix          *readmat = NULL;
1419     real              *eigenvectors;
1420
1421     int                isize = 0, ifsize = 0, iosize = 0;
1422     atom_id           *index = NULL, *fitidx, *outidx;
1423     char              *grpname;
1424     real               rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1425     char               buf[STRLEN], buf1[80], title[STRLEN];
1426     gmx_bool           bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1427
1428     int                method, ncluster = 0;
1429     static const char *methodname[] = {
1430         NULL, "linkage", "jarvis-patrick", "monte-carlo",
1431         "diagonalization", "gromos", NULL
1432     };
1433     enum {
1434         m_null, m_linkage, m_jarvis_patrick,
1435         m_monte_carlo, m_diagonalize, m_gromos, m_nr
1436     };
1437     /* Set colors for plotting: white = zero RMS, black = maximum */
1438     static t_rgb rlo_top  = { 1.0, 1.0, 1.0 };
1439     static t_rgb rhi_top  = { 0.0, 0.0, 0.0 };
1440     static t_rgb rlo_bot  = { 1.0, 1.0, 1.0 };
1441     static t_rgb rhi_bot  = { 0.0, 0.0, 1.0 };
1442     static int   nlevels  = 40, skip = 1;
1443     static real  scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1444     gmx_bool     bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1445     static int   niter    = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1446     static real  kT       = 1e-3;
1447     static int   M        = 10, P = 3;
1448     output_env_t oenv;
1449     gmx_rmpbc_t  gpbc = NULL;
1450
1451     t_pargs      pa[] = {
1452         { "-dista", FALSE, etBOOL, {&bRMSdist},
1453           "Use RMSD of distances instead of RMS deviation" },
1454         { "-nlevels", FALSE, etINT,  {&nlevels},
1455           "Discretize RMSD matrix in this number of levels" },
1456         { "-cutoff", FALSE, etREAL, {&rmsdcut},
1457           "RMSD cut-off (nm) for two structures to be neighbor" },
1458         { "-fit",   FALSE, etBOOL, {&bFit},
1459           "Use least squares fitting before RMSD calculation" },
1460         { "-max",   FALSE, etREAL, {&scalemax},
1461           "Maximum level in RMSD matrix" },
1462         { "-skip",  FALSE, etINT,  {&skip},
1463           "Only analyze every nr-th frame" },
1464         { "-av",    FALSE, etBOOL, {&bAverage},
1465           "Write average iso middle structure for each cluster" },
1466         { "-wcl",   FALSE, etINT,  {&write_ncl},
1467           "Write the structures for this number of clusters to numbered files" },
1468         { "-nst",   FALSE, etINT,  {&write_nst},
1469           "Only write all structures if more than this number of structures per cluster" },
1470         { "-rmsmin", FALSE, etREAL, {&rmsmin},
1471           "minimum rms difference with rest of cluster for writing structures" },
1472         { "-method", FALSE, etENUM, {methodname},
1473           "Method for cluster determination" },
1474         { "-minstruct", FALSE, etINT, {&minstruct},
1475           "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1476         { "-binary", FALSE, etBOOL, {&bBinary},
1477           "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1478           "is given by [TT]-cutoff[tt]" },
1479         { "-M",     FALSE, etINT,  {&M},
1480           "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1481           "0 is use cutoff" },
1482         { "-P",     FALSE, etINT,  {&P},
1483           "Number of identical nearest neighbors required to form a cluster" },
1484         { "-seed",  FALSE, etINT,  {&seed},
1485           "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1486         { "-niter", FALSE, etINT,  {&niter},
1487           "Number of iterations for MC" },
1488         { "-nrandom", FALSE, etINT,  {&nrandom},
1489           "The first iterations for MC may be done complete random, to shuffle the frames" },
1490         { "-kT",    FALSE, etREAL, {&kT},
1491           "Boltzmann weighting factor for Monte Carlo optimization "
1492           "(zero turns off uphill steps)" },
1493         { "-pbc", FALSE, etBOOL,
1494           { &bPBC }, "PBC check" }
1495     };
1496     t_filenm     fnm[] = {
1497         { efTRX, "-f",     NULL,        ffOPTRD },
1498         { efTPS, "-s",     NULL,        ffOPTRD },
1499         { efNDX, NULL,     NULL,        ffOPTRD },
1500         { efXPM, "-dm",   "rmsd",       ffOPTRD },
1501         { efXPM, "-om",   "rmsd-raw",   ffWRITE },
1502         { efXPM, "-o",    "rmsd-clust", ffWRITE },
1503         { efLOG, "-g",    "cluster",    ffWRITE },
1504         { efXVG, "-dist", "rmsd-dist",  ffOPTWR },
1505         { efXVG, "-ev",   "rmsd-eig",   ffOPTWR },
1506         { efXVG, "-conv", "mc-conv",    ffOPTWR },
1507         { efXVG, "-sz",   "clust-size", ffOPTWR},
1508         { efXPM, "-tr",   "clust-trans", ffOPTWR},
1509         { efXVG, "-ntr",  "clust-trans", ffOPTWR},
1510         { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1511         { efTRX, "-cl",   "clusters.pdb", ffOPTWR }
1512     };
1513 #define NFILE asize(fnm)
1514
1515     if (!parse_common_args(&argc, argv,
1516                            PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1517                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1518                            &oenv))
1519     {
1520         return 0;
1521     }
1522
1523     /* parse options */
1524     bReadMat   = opt2bSet("-dm", NFILE, fnm);
1525     bReadTraj  = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1526     if (opt2parg_bSet("-av", asize(pa), pa) ||
1527         opt2parg_bSet("-wcl", asize(pa), pa) ||
1528         opt2parg_bSet("-nst", asize(pa), pa) ||
1529         opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1530         opt2bSet("-cl", NFILE, fnm) )
1531     {
1532         trx_out_fn = opt2fn("-cl", NFILE, fnm);
1533     }
1534     else
1535     {
1536         trx_out_fn = NULL;
1537     }
1538     if (bReadMat && output_env_get_time_factor(oenv) != 1)
1539     {
1540         fprintf(stderr,
1541                 "\nWarning: assuming the time unit in %s is %s\n",
1542                 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1543     }
1544     if (trx_out_fn && !bReadTraj)
1545     {
1546         fprintf(stderr, "\nWarning: "
1547                 "cannot write cluster structures without reading trajectory\n"
1548                 "         ignoring option -cl %s\n", trx_out_fn);
1549     }
1550
1551     method = 1;
1552     while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1553     {
1554         method++;
1555     }
1556     if (method == m_nr)
1557     {
1558         gmx_fatal(FARGS, "Invalid method");
1559     }
1560
1561     bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1562                 method == m_gromos );
1563
1564     /* Open log file */
1565     log = ftp2FILE(efLOG, NFILE, fnm, "w");
1566
1567     fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1568     fprintf(log, "Using %s method for clustering\n", methodname[0]);
1569
1570     /* check input and write parameters to log file */
1571     bUseRmsdCut = FALSE;
1572     if (method == m_jarvis_patrick)
1573     {
1574         bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1575         if ((M < 0) || (M == 1))
1576         {
1577             gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1578         }
1579         if (M < 2)
1580         {
1581             sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1582             bUseRmsdCut = TRUE;
1583         }
1584         else
1585         {
1586             if (P >= M)
1587             {
1588                 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1589             }
1590             if (bJP_RMSD)
1591             {
1592                 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1593                 bUseRmsdCut = TRUE;
1594             }
1595             else
1596             {
1597                 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1598             }
1599         }
1600         ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1601     }
1602     else /* method != m_jarvis */
1603     {
1604         bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1605     }
1606     if (bUseRmsdCut && method != m_jarvis_patrick)
1607     {
1608         fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1609     }
1610     if (method == m_monte_carlo)
1611     {
1612         fprintf(log, "Using %d iterations\n", niter);
1613     }
1614
1615     if (skip < 1)
1616     {
1617         gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1618     }
1619
1620     /* get input */
1621     if (bReadTraj)
1622     {
1623         /* don't read mass-database as masses (and top) are not used */
1624         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1625                       TRUE);
1626         if (bPBC)
1627         {
1628             gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1629         }
1630
1631         fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1632                 bReadMat ? "" : " and RMSD calculation");
1633         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1634                   1, &ifsize, &fitidx, &grpname);
1635         if (trx_out_fn)
1636         {
1637             fprintf(stderr, "\nSelect group for output:\n");
1638             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1639                       1, &iosize, &outidx, &grpname);
1640             /* merge and convert both index groups: */
1641             /* first copy outidx to index. let outidx refer to elements in index */
1642             snew(index, iosize);
1643             isize = iosize;
1644             for (i = 0; i < iosize; i++)
1645             {
1646                 index[i]  = outidx[i];
1647                 outidx[i] = i;
1648             }
1649             /* now lookup elements from fitidx in index, add them if necessary
1650                and also let fitidx refer to elements in index */
1651             for (i = 0; i < ifsize; i++)
1652             {
1653                 j = 0;
1654                 while (j < isize && index[j] != fitidx[i])
1655                 {
1656                     j++;
1657                 }
1658                 if (j >= isize)
1659                 {
1660                     /* slow this way, but doesn't matter much */
1661                     isize++;
1662                     srenew(index, isize);
1663                 }
1664                 index[j]  = fitidx[i];
1665                 fitidx[i] = j;
1666             }
1667         }
1668         else /* !trx_out_fn */
1669         {
1670             isize = ifsize;
1671             snew(index, isize);
1672             for (i = 0; i < ifsize; i++)
1673             {
1674                 index[i]  = fitidx[i];
1675                 fitidx[i] = i;
1676             }
1677         }
1678     }
1679
1680     if (bReadTraj)
1681     {
1682         /* Loop over first coordinate file */
1683         fn = opt2fn("-f", NFILE, fnm);
1684
1685         xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1686         output_env_conv_times(oenv, nf, time);
1687         if (!bRMSdist || bAnalyze)
1688         {
1689             /* Center all frames on zero */
1690             snew(mass, isize);
1691             for (i = 0; i < ifsize; i++)
1692             {
1693                 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1694             }
1695             if (bFit)
1696             {
1697                 for (i = 0; i < nf; i++)
1698                 {
1699                     reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1700                 }
1701             }
1702         }
1703         if (bPBC)
1704         {
1705             gmx_rmpbc_done(gpbc);
1706         }
1707     }
1708
1709     if (bReadMat)
1710     {
1711         fprintf(stderr, "Reading rms distance matrix ");
1712         read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1713         fprintf(stderr, "\n");
1714         if (readmat[0].nx != readmat[0].ny)
1715         {
1716             gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1717                       readmat[0].nx, readmat[0].ny);
1718         }
1719         if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1720         {
1721             gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1722                       "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1723         }
1724
1725         nf = readmat[0].nx;
1726         sfree(time);
1727         time        = readmat[0].axis_x;
1728         time_invfac = output_env_get_time_invfactor(oenv);
1729         for (i = 0; i < nf; i++)
1730         {
1731             time[i] *= time_invfac;
1732         }
1733
1734         rms = init_mat(readmat[0].nx, method == m_diagonalize);
1735         convert_mat(&(readmat[0]), rms);
1736
1737         nlevels = readmat[0].nmap;
1738     }
1739     else   /* !bReadMat */
1740     {
1741         rms  = init_mat(nf, method == m_diagonalize);
1742         nrms = ((gmx_large_int_t)nf*((gmx_large_int_t)nf-1))/2;
1743         if (!bRMSdist)
1744         {
1745             fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1746             /* Initialize work array */
1747             snew(x1, isize);
1748             for (i1 = 0; i1 < nf; i1++)
1749             {
1750                 for (i2 = i1+1; i2 < nf; i2++)
1751                 {
1752                     for (i = 0; i < isize; i++)
1753                     {
1754                         copy_rvec(xx[i1][i], x1[i]);
1755                     }
1756                     if (bFit)
1757                     {
1758                         do_fit(isize, mass, xx[i2], x1);
1759                     }
1760                     rmsd = rmsdev(isize, mass, xx[i2], x1);
1761                     set_mat_entry(rms, i1, i2, rmsd);
1762                 }
1763                 nrms -= (gmx_large_int_t) (nf-i1-1);
1764                 fprintf(stderr, "\r# RMSD calculations left: "gmx_large_int_pfmt"   ", nrms);
1765             }
1766             sfree(x1);
1767         }
1768         else /* bRMSdist */
1769         {
1770             fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1771
1772             /* Initiate work arrays */
1773             snew(d1, isize);
1774             snew(d2, isize);
1775             for (i = 0; (i < isize); i++)
1776             {
1777                 snew(d1[i], isize);
1778                 snew(d2[i], isize);
1779             }
1780             for (i1 = 0; i1 < nf ; i1++)
1781             {
1782                 calc_dist(isize, xx[i1], d1);
1783                 for (i2 = i1+1; (i2 < nf); i2++)
1784                 {
1785                     calc_dist(isize, xx[i2], d2);
1786                     set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1787                 }
1788                 nrms -= (nf-i1-1);
1789                 fprintf(stderr, "\r# RMSD calculations left: "gmx_large_int_pfmt"   ", nrms);
1790             }
1791             /* Clean up work arrays */
1792             for (i = 0; (i < isize); i++)
1793             {
1794                 sfree(d1[i]);
1795                 sfree(d2[i]);
1796             }
1797             sfree(d1);
1798             sfree(d2);
1799         }
1800         fprintf(stderr, "\n\n");
1801     }
1802     ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1803                 rms->minrms, rms->maxrms);
1804     ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1805     ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1806     ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1807     if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1808     {
1809         fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1810                 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1811     }
1812     if (bAnalyze && (rmsmin < rms->minrms) )
1813     {
1814         fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1815                 rmsmin, rms->minrms);
1816     }
1817     if (bAnalyze && (rmsmin > rmsdcut) )
1818     {
1819         fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1820                 rmsmin, rmsdcut);
1821     }
1822
1823     /* Plot the rmsd distribution */
1824     rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1825
1826     if (bBinary)
1827     {
1828         for (i1 = 0; (i1 < nf); i1++)
1829         {
1830             for (i2 = 0; (i2 < nf); i2++)
1831             {
1832                 if (rms->mat[i1][i2] < rmsdcut)
1833                 {
1834                     rms->mat[i1][i2] = 0;
1835                 }
1836                 else
1837                 {
1838                     rms->mat[i1][i2] = 1;
1839                 }
1840             }
1841         }
1842     }
1843
1844     snew(clust.cl, nf);
1845     switch (method)
1846     {
1847         case m_linkage:
1848             /* Now sort the matrix and write it out again */
1849             gather(rms, rmsdcut, &clust);
1850             break;
1851         case m_diagonalize:
1852             /* Do a diagonalization */
1853             snew(eigenvalues, nf);
1854             snew(eigenvectors, nf*nf);
1855             memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1856             eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1857             sfree(eigenvectors);
1858
1859             fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1860                           "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1861             for (i = 0; (i < nf); i++)
1862             {
1863                 fprintf(fp, "%10d  %10g\n", i, eigenvalues[i]);
1864             }
1865             ffclose(fp);
1866             break;
1867         case m_monte_carlo:
1868             orig     = init_mat(rms->nn, FALSE);
1869             orig->nn = rms->nn;
1870             copy_t_mat(orig, rms);
1871             mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1872                         opt2fn_null("-conv", NFILE, fnm), oenv);
1873             break;
1874         case m_jarvis_patrick:
1875             jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1876             break;
1877         case m_gromos:
1878             gromos(rms->nn, rms->mat, rmsdcut, &clust);
1879             break;
1880         default:
1881             gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1882     }
1883
1884     if (method == m_monte_carlo || method == m_diagonalize)
1885     {
1886         fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1887                 mat_energy(rms));
1888     }
1889
1890     if (bAnalyze)
1891     {
1892         if (minstruct > 1)
1893         {
1894             ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1895         }
1896         else
1897         {
1898             mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1899         }
1900         init_t_atoms(&useatoms, isize, FALSE);
1901         snew(usextps, isize);
1902         useatoms.resinfo = top.atoms.resinfo;
1903         for (i = 0; i < isize; i++)
1904         {
1905             useatoms.atomname[i]    = top.atoms.atomname[index[i]];
1906             useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1907             useatoms.nres           = max(useatoms.nres, useatoms.atom[i].resind+1);
1908             copy_rvec(xtps[index[i]], usextps[i]);
1909         }
1910         useatoms.nr = isize;
1911         analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1912                          ifsize, fitidx, iosize, outidx,
1913                          bReadTraj ? trx_out_fn : NULL,
1914                          opt2fn_null("-sz", NFILE, fnm),
1915                          opt2fn_null("-tr", NFILE, fnm),
1916                          opt2fn_null("-ntr", NFILE, fnm),
1917                          opt2fn_null("-clid", NFILE, fnm),
1918                          bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1919                          rlo_bot, rhi_bot, oenv);
1920     }
1921     ffclose(log);
1922
1923     if (bBinary && !bAnalyze)
1924     {
1925         /* Make the clustering visible */
1926         for (i2 = 0; (i2 < nf); i2++)
1927         {
1928             for (i1 = i2+1; (i1 < nf); i1++)
1929             {
1930                 if (rms->mat[i1][i2])
1931                 {
1932                     rms->mat[i1][i2] = rms->maxrms;
1933                 }
1934             }
1935         }
1936     }
1937
1938     fp = opt2FILE("-o", NFILE, fnm, "w");
1939     fprintf(stderr, "Writing rms distance/clustering matrix ");
1940     if (bReadMat)
1941     {
1942         write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1943                   readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1944                   rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1945     }
1946     else
1947     {
1948         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1949         sprintf(title, "RMS%sDeviation / Cluster Index",
1950                 bRMSdist ? " Distance " : " ");
1951         if (minstruct > 1)
1952         {
1953             write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1954                             nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1955                             rlo_top, rhi_top, 0.0, (real) ncluster,
1956                             &ncluster, TRUE, rlo_bot, rhi_bot);
1957         }
1958         else
1959         {
1960             write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1961                       nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1962                       rlo_top, rhi_top, &nlevels);
1963         }
1964     }
1965     fprintf(stderr, "\n");
1966     ffclose(fp);
1967     if (NULL != orig)
1968     {
1969         fp = opt2FILE("-om", NFILE, fnm, "w");
1970         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1971         sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1972         write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1973                   nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1974                   rlo_top, rhi_top, &nlevels);
1975         ffclose(fp);
1976         done_mat(&orig);
1977         sfree(orig);
1978     }
1979     /* now show what we've done */
1980     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1981     do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1982     if (method == m_diagonalize)
1983     {
1984         do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1985     }
1986     do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1987     if (bAnalyze)
1988     {
1989         do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1990         do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1991         do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1992     }
1993     do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);
1994
1995     return 0;
1996 }