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