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