src/gromacs/gmxana/gmx_cluster.cpp

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