src/gromacs/gmxana/gmx_cluster.c

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