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