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