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