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                i, i1, i2, j, nf, nrms;
1407
1408     matrix             box;
1409     rvec              *xtps, *usextps, *x1, **xx = NULL;
1410     const char        *fn, *trx_out_fn;
1411     t_clusters         clust;
1412     t_mat             *rms, *orig = NULL;
1413     real              *eigenvalues;
1414     t_topology         top;
1415     int                ePBC;
1416     t_atoms            useatoms;
1417     t_matrix          *readmat = NULL;
1418     real              *eigenvectors;
1419
1420     int                isize = 0, ifsize = 0, iosize = 0;
1421     atom_id           *index = NULL, *fitidx, *outidx;
1422     char              *grpname;
1423     real               rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1424     char               buf[STRLEN], buf1[80], title[STRLEN];
1425     gmx_bool           bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1426
1427     int                method, ncluster = 0;
1428     static const char *methodname[] = {
1429         NULL, "linkage", "jarvis-patrick", "monte-carlo",
1430         "diagonalization", "gromos", NULL
1431     };
1432     enum {
1433         m_null, m_linkage, m_jarvis_patrick,
1434         m_monte_carlo, m_diagonalize, m_gromos, m_nr
1435     };
1436     /* Set colors for plotting: white = zero RMS, black = maximum */
1437     static t_rgb rlo_top  = { 1.0, 1.0, 1.0 };
1438     static t_rgb rhi_top  = { 0.0, 0.0, 0.0 };
1439     static t_rgb rlo_bot  = { 1.0, 1.0, 1.0 };
1440     static t_rgb rhi_bot  = { 0.0, 0.0, 1.0 };
1441     static int   nlevels  = 40, skip = 1;
1442     static real  scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1443     gmx_bool     bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1444     static int   niter    = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1445     static real  kT       = 1e-3;
1446     static int   M        = 10, P = 3;
1447     output_env_t oenv;
1448     gmx_rmpbc_t  gpbc = NULL;
1449
1450     t_pargs      pa[] = {
1451         { "-dista", FALSE, etBOOL, {&bRMSdist},
1452           "Use RMSD of distances instead of RMS deviation" },
1453         { "-nlevels", FALSE, etINT,  {&nlevels},
1454           "Discretize RMSD matrix in this number of levels" },
1455         { "-cutoff", FALSE, etREAL, {&rmsdcut},
1456           "RMSD cut-off (nm) for two structures to be neighbor" },
1457         { "-fit",   FALSE, etBOOL, {&bFit},
1458           "Use least squares fitting before RMSD calculation" },
1459         { "-max",   FALSE, etREAL, {&scalemax},
1460           "Maximum level in RMSD matrix" },
1461         { "-skip",  FALSE, etINT,  {&skip},
1462           "Only analyze every nr-th frame" },
1463         { "-av",    FALSE, etBOOL, {&bAverage},
1464           "Write average iso middle structure for each cluster" },
1465         { "-wcl",   FALSE, etINT,  {&write_ncl},
1466           "Write the structures for this number of clusters to numbered files" },
1467         { "-nst",   FALSE, etINT,  {&write_nst},
1468           "Only write all structures if more than this number of structures per cluster" },
1469         { "-rmsmin", FALSE, etREAL, {&rmsmin},
1470           "minimum rms difference with rest of cluster for writing structures" },
1471         { "-method", FALSE, etENUM, {methodname},
1472           "Method for cluster determination" },
1473         { "-minstruct", FALSE, etINT, {&minstruct},
1474           "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1475         { "-binary", FALSE, etBOOL, {&bBinary},
1476           "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1477           "is given by [TT]-cutoff[tt]" },
1478         { "-M",     FALSE, etINT,  {&M},
1479           "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1480           "0 is use cutoff" },
1481         { "-P",     FALSE, etINT,  {&P},
1482           "Number of identical nearest neighbors required to form a cluster" },
1483         { "-seed",  FALSE, etINT,  {&seed},
1484           "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1485         { "-niter", FALSE, etINT,  {&niter},
1486           "Number of iterations for MC" },
1487         { "-nrandom", FALSE, etINT,  {&nrandom},
1488           "The first iterations for MC may be done complete random, to shuffle the frames" },
1489         { "-kT",    FALSE, etREAL, {&kT},
1490           "Boltzmann weighting factor for Monte Carlo optimization "
1491           "(zero turns off uphill steps)" },
1492         { "-pbc", FALSE, etBOOL,
1493           { &bPBC }, "PBC check" }
1494     };
1495     t_filenm     fnm[] = {
1496         { efTRX, "-f",     NULL,        ffOPTRD },
1497         { efTPS, "-s",     NULL,        ffOPTRD },
1498         { efNDX, NULL,     NULL,        ffOPTRD },
1499         { efXPM, "-dm",   "rmsd",       ffOPTRD },
1500         { efXPM, "-om",   "rmsd-raw",   ffWRITE },
1501         { efXPM, "-o",    "rmsd-clust", ffWRITE },
1502         { efLOG, "-g",    "cluster",    ffWRITE },
1503         { efXVG, "-dist", "rmsd-dist",  ffOPTWR },
1504         { efXVG, "-ev",   "rmsd-eig",   ffOPTWR },
1505         { efXVG, "-conv", "mc-conv",    ffOPTWR },
1506         { efXVG, "-sz",   "clust-size", ffOPTWR},
1507         { efXPM, "-tr",   "clust-trans", ffOPTWR},
1508         { efXVG, "-ntr",  "clust-trans", ffOPTWR},
1509         { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1510         { efTRX, "-cl",   "clusters.pdb", ffOPTWR }
1511     };
1512 #define NFILE asize(fnm)
1513
1514     if (!parse_common_args(&argc, argv,
1515                            PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1516                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1517                            &oenv))
1518     {
1519         return 0;
1520     }
1521
1522     /* parse options */
1523     bReadMat   = opt2bSet("-dm", NFILE, fnm);
1524     bReadTraj  = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1525     if (opt2parg_bSet("-av", asize(pa), pa) ||
1526         opt2parg_bSet("-wcl", asize(pa), pa) ||
1527         opt2parg_bSet("-nst", asize(pa), pa) ||
1528         opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1529         opt2bSet("-cl", NFILE, fnm) )
1530     {
1531         trx_out_fn = opt2fn("-cl", NFILE, fnm);
1532     }
1533     else
1534     {
1535         trx_out_fn = NULL;
1536     }
1537     if (bReadMat && output_env_get_time_factor(oenv) != 1)
1538     {
1539         fprintf(stderr,
1540                 "\nWarning: assuming the time unit in %s is %s\n",
1541                 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1542     }
1543     if (trx_out_fn && !bReadTraj)
1544     {
1545         fprintf(stderr, "\nWarning: "
1546                 "cannot write cluster structures without reading trajectory\n"
1547                 "         ignoring option -cl %s\n", trx_out_fn);
1548     }
1549
1550     method = 1;
1551     while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1552     {
1553         method++;
1554     }
1555     if (method == m_nr)
1556     {
1557         gmx_fatal(FARGS, "Invalid method");
1558     }
1559
1560     bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1561                 method == m_gromos );
1562
1563     /* Open log file */
1564     log = ftp2FILE(efLOG, NFILE, fnm, "w");
1565
1566     fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1567     fprintf(log, "Using %s method for clustering\n", methodname[0]);
1568
1569     /* check input and write parameters to log file */
1570     bUseRmsdCut = FALSE;
1571     if (method == m_jarvis_patrick)
1572     {
1573         bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1574         if ((M < 0) || (M == 1))
1575         {
1576             gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1577         }
1578         if (M < 2)
1579         {
1580             sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1581             bUseRmsdCut = TRUE;
1582         }
1583         else
1584         {
1585             if (P >= M)
1586             {
1587                 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1588             }
1589             if (bJP_RMSD)
1590             {
1591                 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1592                 bUseRmsdCut = TRUE;
1593             }
1594             else
1595             {
1596                 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1597             }
1598         }
1599         ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1600     }
1601     else /* method != m_jarvis */
1602     {
1603         bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1604     }
1605     if (bUseRmsdCut && method != m_jarvis_patrick)
1606     {
1607         fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1608     }
1609     if (method == m_monte_carlo)
1610     {
1611         fprintf(log, "Using %d iterations\n", niter);
1612     }
1613
1614     if (skip < 1)
1615     {
1616         gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1617     }
1618
1619     /* get input */
1620     if (bReadTraj)
1621     {
1622         /* don't read mass-database as masses (and top) are not used */
1623         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1624                       TRUE);
1625         if (bPBC)
1626         {
1627             gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1628         }
1629
1630         fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1631                 bReadMat ? "" : " and RMSD calculation");
1632         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1633                   1, &ifsize, &fitidx, &grpname);
1634         if (trx_out_fn)
1635         {
1636             fprintf(stderr, "\nSelect group for output:\n");
1637             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1638                       1, &iosize, &outidx, &grpname);
1639             /* merge and convert both index groups: */
1640             /* first copy outidx to index. let outidx refer to elements in index */
1641             snew(index, iosize);
1642             isize = iosize;
1643             for (i = 0; i < iosize; i++)
1644             {
1645                 index[i]  = outidx[i];
1646                 outidx[i] = i;
1647             }
1648             /* now lookup elements from fitidx in index, add them if necessary
1649                and also let fitidx refer to elements in index */
1650             for (i = 0; i < ifsize; i++)
1651             {
1652                 j = 0;
1653                 while (j < isize && index[j] != fitidx[i])
1654                 {
1655                     j++;
1656                 }
1657                 if (j >= isize)
1658                 {
1659                     /* slow this way, but doesn't matter much */
1660                     isize++;
1661                     srenew(index, isize);
1662                 }
1663                 index[j]  = fitidx[i];
1664                 fitidx[i] = j;
1665             }
1666         }
1667         else /* !trx_out_fn */
1668         {
1669             isize = ifsize;
1670             snew(index, isize);
1671             for (i = 0; i < ifsize; i++)
1672             {
1673                 index[i]  = fitidx[i];
1674                 fitidx[i] = i;
1675             }
1676         }
1677     }
1678     /* Initiate arrays */
1679     snew(d1, isize);
1680     snew(d2, isize);
1681     for (i = 0; (i < isize); i++)
1682     {
1683         snew(d1[i], isize);
1684         snew(d2[i], isize);
1685     }
1686
1687     if (bReadTraj)
1688     {
1689         /* Loop over first coordinate file */
1690         fn = opt2fn("-f", NFILE, fnm);
1691
1692         xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1693         output_env_conv_times(oenv, nf, time);
1694         if (!bRMSdist || bAnalyze)
1695         {
1696             /* Center all frames on zero */
1697             snew(mass, isize);
1698             for (i = 0; i < ifsize; i++)
1699             {
1700                 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1701             }
1702             if (bFit)
1703             {
1704                 for (i = 0; i < nf; i++)
1705                 {
1706                     reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1707                 }
1708             }
1709         }
1710         if (bPBC)
1711         {
1712             gmx_rmpbc_done(gpbc);
1713         }
1714     }
1715
1716     if (bReadMat)
1717     {
1718         fprintf(stderr, "Reading rms distance matrix ");
1719         read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1720         fprintf(stderr, "\n");
1721         if (readmat[0].nx != readmat[0].ny)
1722         {
1723             gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1724                       readmat[0].nx, readmat[0].ny);
1725         }
1726         if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1727         {
1728             gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1729                       "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1730         }
1731
1732         nf = readmat[0].nx;
1733         sfree(time);
1734         time        = readmat[0].axis_x;
1735         time_invfac = output_env_get_time_invfactor(oenv);
1736         for (i = 0; i < nf; i++)
1737         {
1738             time[i] *= time_invfac;
1739         }
1740
1741         rms = init_mat(readmat[0].nx, method == m_diagonalize);
1742         convert_mat(&(readmat[0]), rms);
1743
1744         nlevels = readmat[0].nmap;
1745     }
1746     else   /* !bReadMat */
1747     {
1748         rms  = init_mat(nf, method == m_diagonalize);
1749         nrms = (nf*(nf-1))/2;
1750         if (!bRMSdist)
1751         {
1752             fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1753             snew(x1, isize);
1754             for (i1 = 0; (i1 < nf); i1++)
1755             {
1756                 for (i2 = i1+1; (i2 < nf); i2++)
1757                 {
1758                     for (i = 0; i < isize; i++)
1759                     {
1760                         copy_rvec(xx[i1][i], x1[i]);
1761                     }
1762                     if (bFit)
1763                     {
1764                         do_fit(isize, mass, xx[i2], x1);
1765                     }
1766                     rmsd = rmsdev(isize, mass, xx[i2], x1);
1767                     set_mat_entry(rms, i1, i2, rmsd);
1768                 }
1769                 nrms -= (nf-i1-1);
1770                 fprintf(stderr, "\r# RMSD calculations left: %d   ", nrms);
1771             }
1772         }
1773         else /* bRMSdist */
1774         {
1775             fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1776             for (i1 = 0; (i1 < nf); i1++)
1777             {
1778                 calc_dist(isize, xx[i1], d1);
1779                 for (i2 = i1+1; (i2 < nf); i2++)
1780                 {
1781                     calc_dist(isize, xx[i2], d2);
1782                     set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1783                 }
1784                 nrms -= (nf-i1-1);
1785                 fprintf(stderr, "\r# RMSD calculations left: %d   ", nrms);
1786             }
1787         }
1788         fprintf(stderr, "\n\n");
1789     }
1790     ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1791                 rms->minrms, rms->maxrms);
1792     ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1793     ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1794     ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1795     if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1796     {
1797         fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1798                 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1799     }
1800     if (bAnalyze && (rmsmin < rms->minrms) )
1801     {
1802         fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1803                 rmsmin, rms->minrms);
1804     }
1805     if (bAnalyze && (rmsmin > rmsdcut) )
1806     {
1807         fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1808                 rmsmin, rmsdcut);
1809     }
1810
1811     /* Plot the rmsd distribution */
1812     rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1813
1814     if (bBinary)
1815     {
1816         for (i1 = 0; (i1 < nf); i1++)
1817         {
1818             for (i2 = 0; (i2 < nf); i2++)
1819             {
1820                 if (rms->mat[i1][i2] < rmsdcut)
1821                 {
1822                     rms->mat[i1][i2] = 0;
1823                 }
1824                 else
1825                 {
1826                     rms->mat[i1][i2] = 1;
1827                 }
1828             }
1829         }
1830     }
1831
1832     snew(clust.cl, nf);
1833     switch (method)
1834     {
1835         case m_linkage:
1836             /* Now sort the matrix and write it out again */
1837             gather(rms, rmsdcut, &clust);
1838             break;
1839         case m_diagonalize:
1840             /* Do a diagonalization */
1841             snew(eigenvalues, nf);
1842             snew(eigenvectors, nf*nf);
1843             memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1844             eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1845             sfree(eigenvectors);
1846
1847             fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1848                           "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1849             for (i = 0; (i < nf); i++)
1850             {
1851                 fprintf(fp, "%10d  %10g\n", i, eigenvalues[i]);
1852             }
1853             ffclose(fp);
1854             break;
1855         case m_monte_carlo:
1856             orig     = init_mat(rms->nn, FALSE);
1857             orig->nn = rms->nn;
1858             copy_t_mat(orig, rms);
1859             mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1860                         opt2fn_null("-conv", NFILE, fnm), oenv);
1861             break;
1862         case m_jarvis_patrick:
1863             jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1864             break;
1865         case m_gromos:
1866             gromos(rms->nn, rms->mat, rmsdcut, &clust);
1867             break;
1868         default:
1869             gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1870     }
1871
1872     if (method == m_monte_carlo || method == m_diagonalize)
1873     {
1874         fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1875                 mat_energy(rms));
1876     }
1877
1878     if (bAnalyze)
1879     {
1880         if (minstruct > 1)
1881         {
1882             ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1883         }
1884         else
1885         {
1886             mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1887         }
1888         init_t_atoms(&useatoms, isize, FALSE);
1889         snew(usextps, isize);
1890         useatoms.resinfo = top.atoms.resinfo;
1891         for (i = 0; i < isize; i++)
1892         {
1893             useatoms.atomname[i]    = top.atoms.atomname[index[i]];
1894             useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1895             useatoms.nres           = max(useatoms.nres, useatoms.atom[i].resind+1);
1896             copy_rvec(xtps[index[i]], usextps[i]);
1897         }
1898         useatoms.nr = isize;
1899         analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1900                          ifsize, fitidx, iosize, outidx,
1901                          bReadTraj ? trx_out_fn : NULL,
1902                          opt2fn_null("-sz", NFILE, fnm),
1903                          opt2fn_null("-tr", NFILE, fnm),
1904                          opt2fn_null("-ntr", NFILE, fnm),
1905                          opt2fn_null("-clid", NFILE, fnm),
1906                          bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1907                          rlo_bot, rhi_bot, oenv);
1908     }
1909     ffclose(log);
1910
1911     if (bBinary && !bAnalyze)
1912     {
1913         /* Make the clustering visible */
1914         for (i2 = 0; (i2 < nf); i2++)
1915         {
1916             for (i1 = i2+1; (i1 < nf); i1++)
1917             {
1918                 if (rms->mat[i1][i2])
1919                 {
1920                     rms->mat[i1][i2] = rms->maxrms;
1921                 }
1922             }
1923         }
1924     }
1925
1926     fp = opt2FILE("-o", NFILE, fnm, "w");
1927     fprintf(stderr, "Writing rms distance/clustering matrix ");
1928     if (bReadMat)
1929     {
1930         write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1931                   readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1932                   rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1933     }
1934     else
1935     {
1936         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1937         sprintf(title, "RMS%sDeviation / Cluster Index",
1938                 bRMSdist ? " Distance " : " ");
1939         if (minstruct > 1)
1940         {
1941             write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1942                             nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1943                             rlo_top, rhi_top, 0.0, (real) ncluster,
1944                             &ncluster, TRUE, rlo_bot, rhi_bot);
1945         }
1946         else
1947         {
1948             write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1949                       nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1950                       rlo_top, rhi_top, &nlevels);
1951         }
1952     }
1953     fprintf(stderr, "\n");
1954     ffclose(fp);
1955     if (NULL != orig)
1956     {
1957         fp = opt2FILE("-om", NFILE, fnm, "w");
1958         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1959         sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1960         write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1961                   nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1962                   rlo_top, rhi_top, &nlevels);
1963         ffclose(fp);
1964         done_mat(&orig);
1965         sfree(orig);
1966     }
1967     /* now show what we've done */
1968     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1969     do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1970     if (method == m_diagonalize)
1971     {
1972         do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1973     }
1974     do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1975     if (bAnalyze)
1976     {
1977         do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1978         do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1979         do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1980     }
1981     do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);
1982
1983     return 0;
1984 }