src/gromacs/gmxana/gmx_cluster.c

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