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     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     /* parse options */
1520     bReadMat   = opt2bSet("-dm", NFILE, fnm);
1521     bReadTraj  = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1522     if (opt2parg_bSet("-av", asize(pa), pa) ||
1523         opt2parg_bSet("-wcl", asize(pa), pa) ||
1524         opt2parg_bSet("-nst", asize(pa), pa) ||
1525         opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1526         opt2bSet("-cl", NFILE, fnm) )
1527     {
1528         trx_out_fn = opt2fn("-cl", NFILE, fnm);
1529     }
1530     else
1531     {
1532         trx_out_fn = NULL;
1533     }
1534     if (bReadMat && output_env_get_time_factor(oenv) != 1)
1535     {
1536         fprintf(stderr,
1537                 "\nWarning: assuming the time unit in %s is %s\n",
1538                 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1539     }
1540     if (trx_out_fn && !bReadTraj)
1541     {
1542         fprintf(stderr, "\nWarning: "
1543                 "cannot write cluster structures without reading trajectory\n"
1544                 "         ignoring option -cl %s\n", trx_out_fn);
1545     }
1546
1547     method = 1;
1548     while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1549     {
1550         method++;
1551     }
1552     if (method == m_nr)
1553     {
1554         gmx_fatal(FARGS, "Invalid method");
1555     }
1556
1557     bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1558                 method == m_gromos );
1559
1560     /* Open log file */
1561     log = ftp2FILE(efLOG, NFILE, fnm, "w");
1562
1563     fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1564     fprintf(log, "Using %s method for clustering\n", methodname[0]);
1565
1566     /* check input and write parameters to log file */
1567     bUseRmsdCut = FALSE;
1568     if (method == m_jarvis_patrick)
1569     {
1570         bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1571         if ((M < 0) || (M == 1))
1572         {
1573             gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1574         }
1575         if (M < 2)
1576         {
1577             sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1578             bUseRmsdCut = TRUE;
1579         }
1580         else
1581         {
1582             if (P >= M)
1583             {
1584                 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1585             }
1586             if (bJP_RMSD)
1587             {
1588                 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1589                 bUseRmsdCut = TRUE;
1590             }
1591             else
1592             {
1593                 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1594             }
1595         }
1596         ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1597     }
1598     else /* method != m_jarvis */
1599     {
1600         bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1601     }
1602     if (bUseRmsdCut && method != m_jarvis_patrick)
1603     {
1604         fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1605     }
1606     if (method == m_monte_carlo)
1607     {
1608         fprintf(log, "Using %d iterations\n", niter);
1609     }
1610
1611     if (skip < 1)
1612     {
1613         gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1614     }
1615
1616     /* get input */
1617     if (bReadTraj)
1618     {
1619         /* don't read mass-database as masses (and top) are not used */
1620         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1621                       TRUE);
1622         if (bPBC)
1623         {
1624             gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1625         }
1626
1627         fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1628                 bReadMat ? "" : " and RMSD calculation");
1629         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1630                   1, &ifsize, &fitidx, &grpname);
1631         if (trx_out_fn)
1632         {
1633             fprintf(stderr, "\nSelect group for output:\n");
1634             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1635                       1, &iosize, &outidx, &grpname);
1636             /* merge and convert both index groups: */
1637             /* first copy outidx to index. let outidx refer to elements in index */
1638             snew(index, iosize);
1639             isize = iosize;
1640             for (i = 0; i < iosize; i++)
1641             {
1642                 index[i]  = outidx[i];
1643                 outidx[i] = i;
1644             }
1645             /* now lookup elements from fitidx in index, add them if necessary
1646                and also let fitidx refer to elements in index */
1647             for (i = 0; i < ifsize; i++)
1648             {
1649                 j = 0;
1650                 while (j < isize && index[j] != fitidx[i])
1651                 {
1652                     j++;
1653                 }
1654                 if (j >= isize)
1655                 {
1656                     /* slow this way, but doesn't matter much */
1657                     isize++;
1658                     srenew(index, isize);
1659                 }
1660                 index[j]  = fitidx[i];
1661                 fitidx[i] = j;
1662             }
1663         }
1664         else /* !trx_out_fn */
1665         {
1666             isize = ifsize;
1667             snew(index, isize);
1668             for (i = 0; i < ifsize; i++)
1669             {
1670                 index[i]  = fitidx[i];
1671                 fitidx[i] = i;
1672             }
1673         }
1674     }
1675     /* Initiate arrays */
1676     snew(d1, isize);
1677     snew(d2, isize);
1678     for (i = 0; (i < isize); i++)
1679     {
1680         snew(d1[i], isize);
1681         snew(d2[i], isize);
1682     }
1683
1684     if (bReadTraj)
1685     {
1686         /* Loop over first coordinate file */
1687         fn = opt2fn("-f", NFILE, fnm);
1688
1689         xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1690         output_env_conv_times(oenv, nf, time);
1691         if (!bRMSdist || bAnalyze)
1692         {
1693             /* Center all frames on zero */
1694             snew(mass, isize);
1695             for (i = 0; i < ifsize; i++)
1696             {
1697                 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1698             }
1699             if (bFit)
1700             {
1701                 for (i = 0; i < nf; i++)
1702                 {
1703                     reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1704                 }
1705             }
1706         }
1707         if (bPBC)
1708         {
1709             gmx_rmpbc_done(gpbc);
1710         }
1711     }
1712
1713     if (bReadMat)
1714     {
1715         fprintf(stderr, "Reading rms distance matrix ");
1716         read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1717         fprintf(stderr, "\n");
1718         if (readmat[0].nx != readmat[0].ny)
1719         {
1720             gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1721                       readmat[0].nx, readmat[0].ny);
1722         }
1723         if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1724         {
1725             gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1726                       "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1727         }
1728
1729         nf = readmat[0].nx;
1730         sfree(time);
1731         time        = readmat[0].axis_x;
1732         time_invfac = output_env_get_time_invfactor(oenv);
1733         for (i = 0; i < nf; i++)
1734         {
1735             time[i] *= time_invfac;
1736         }
1737
1738         rms = init_mat(readmat[0].nx, method == m_diagonalize);
1739         convert_mat(&(readmat[0]), rms);
1740
1741         nlevels = readmat[0].nmap;
1742     }
1743     else   /* !bReadMat */
1744     {
1745         rms  = init_mat(nf, method == m_diagonalize);
1746         nrms = (nf*(nf-1))/2;
1747         if (!bRMSdist)
1748         {
1749             fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1750             snew(x1, isize);
1751             for (i1 = 0; (i1 < nf); i1++)
1752             {
1753                 for (i2 = i1+1; (i2 < nf); i2++)
1754                 {
1755                     for (i = 0; i < isize; i++)
1756                     {
1757                         copy_rvec(xx[i1][i], x1[i]);
1758                     }
1759                     if (bFit)
1760                     {
1761                         do_fit(isize, mass, xx[i2], x1);
1762                     }
1763                     rmsd = rmsdev(isize, mass, xx[i2], x1);
1764                     set_mat_entry(rms, i1, i2, rmsd);
1765                 }
1766                 nrms -= (nf-i1-1);
1767                 fprintf(stderr, "\r# RMSD calculations left: %d   ", nrms);
1768             }
1769         }
1770         else /* bRMSdist */
1771         {
1772             fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1773             for (i1 = 0; (i1 < nf); i1++)
1774             {
1775                 calc_dist(isize, xx[i1], d1);
1776                 for (i2 = i1+1; (i2 < nf); i2++)
1777                 {
1778                     calc_dist(isize, xx[i2], d2);
1779                     set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1780                 }
1781                 nrms -= (nf-i1-1);
1782                 fprintf(stderr, "\r# RMSD calculations left: %d   ", nrms);
1783             }
1784         }
1785         fprintf(stderr, "\n\n");
1786     }
1787     ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1788                 rms->minrms, rms->maxrms);
1789     ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1790     ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1791     ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1792     if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1793     {
1794         fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1795                 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1796     }
1797     if (bAnalyze && (rmsmin < rms->minrms) )
1798     {
1799         fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1800                 rmsmin, rms->minrms);
1801     }
1802     if (bAnalyze && (rmsmin > rmsdcut) )
1803     {
1804         fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1805                 rmsmin, rmsdcut);
1806     }
1807
1808     /* Plot the rmsd distribution */
1809     rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1810
1811     if (bBinary)
1812     {
1813         for (i1 = 0; (i1 < nf); i1++)
1814         {
1815             for (i2 = 0; (i2 < nf); i2++)
1816             {
1817                 if (rms->mat[i1][i2] < rmsdcut)
1818                 {
1819                     rms->mat[i1][i2] = 0;
1820                 }
1821                 else
1822                 {
1823                     rms->mat[i1][i2] = 1;
1824                 }
1825             }
1826         }
1827     }
1828
1829     snew(clust.cl, nf);
1830     switch (method)
1831     {
1832         case m_linkage:
1833             /* Now sort the matrix and write it out again */
1834             gather(rms, rmsdcut, &clust);
1835             break;
1836         case m_diagonalize:
1837             /* Do a diagonalization */
1838             snew(eigenvalues, nf);
1839             snew(eigenvectors, nf*nf);
1840             memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1841             eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1842             sfree(eigenvectors);
1843
1844             fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1845                           "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1846             for (i = 0; (i < nf); i++)
1847             {
1848                 fprintf(fp, "%10d  %10g\n", i, eigenvalues[i]);
1849             }
1850             ffclose(fp);
1851             break;
1852         case m_monte_carlo:
1853             orig = init_mat(rms->nn, FALSE);
1854             orig->nn = rms->nn;
1855             copy_t_mat(orig, rms);
1856             mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1857                         opt2fn_null("-conv",NFILE,fnm), oenv);
1858             break;
1859         case m_jarvis_patrick:
1860             jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1861             break;
1862         case m_gromos:
1863             gromos(rms->nn, rms->mat, rmsdcut, &clust);
1864             break;
1865         default:
1866             gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1867     }
1868
1869     if (method == m_monte_carlo || method == m_diagonalize)
1870     {
1871         fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1872                 mat_energy(rms));
1873     }
1874
1875     if (bAnalyze)
1876     {
1877         if (minstruct > 1)
1878         {
1879             ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1880         }
1881         else
1882         {
1883             mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1884         }
1885         init_t_atoms(&useatoms, isize, FALSE);
1886         snew(usextps, isize);
1887         useatoms.resinfo = top.atoms.resinfo;
1888         for (i = 0; i < isize; i++)
1889         {
1890             useatoms.atomname[i]    = top.atoms.atomname[index[i]];
1891             useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1892             useatoms.nres           = max(useatoms.nres, useatoms.atom[i].resind+1);
1893             copy_rvec(xtps[index[i]], usextps[i]);
1894         }
1895         useatoms.nr = isize;
1896         analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1897                          ifsize, fitidx, iosize, outidx,
1898                          bReadTraj ? trx_out_fn : NULL,
1899                          opt2fn_null("-sz", NFILE, fnm),
1900                          opt2fn_null("-tr", NFILE, fnm),
1901                          opt2fn_null("-ntr", NFILE, fnm),
1902                          opt2fn_null("-clid", NFILE, fnm),
1903                          bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1904                          rlo_bot, rhi_bot, oenv);
1905     }
1906     ffclose(log);
1907
1908     if (bBinary && !bAnalyze)
1909     {
1910         /* Make the clustering visible */
1911         for (i2 = 0; (i2 < nf); i2++)
1912         {
1913             for (i1 = i2+1; (i1 < nf); i1++)
1914             {
1915                 if (rms->mat[i1][i2])
1916                 {
1917                     rms->mat[i1][i2] = rms->maxrms;
1918                 }
1919             }
1920         }
1921     }
1922
1923     fp = opt2FILE("-o", NFILE, fnm, "w");
1924     fprintf(stderr, "Writing rms distance/clustering matrix ");
1925     if (bReadMat)
1926     {
1927         write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1928                   readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1929                   rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1930     }
1931     else
1932     {
1933         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1934         sprintf(title, "RMS%sDeviation / Cluster Index",
1935                 bRMSdist ? " Distance " : " ");
1936         if (minstruct > 1)
1937         {
1938             write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1939                             nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1940                             rlo_top, rhi_top, 0.0, (real) ncluster,
1941                             &ncluster, TRUE, rlo_bot, rhi_bot);
1942         }
1943         else
1944         {
1945             write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1946                       nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1947                       rlo_top, rhi_top, &nlevels);
1948         }
1949     }
1950     fprintf(stderr, "\n");
1951     ffclose(fp);
1952     if (NULL != orig)
1953     {
1954         fp = opt2FILE("-om", NFILE, fnm, "w");
1955         sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1956         sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1957         write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1958                   nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1959                   rlo_top, rhi_top, &nlevels);
1960         ffclose(fp);
1961         done_mat(&orig);
1962         sfree(orig);
1963     }
1964     /* now show what we've done */
1965     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1966     do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1967     if (method == m_diagonalize)
1968     {
1969         do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1970     }
1971     do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1972     if (bAnalyze)
1973     {
1974         do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1975         do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1976         do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1977     }
1978     do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);
1979
1980     return 0;
1981 }