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,2018,2019, 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 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 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 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 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 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 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 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 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 /= gmx::exactDiv(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] = static_cast<real>(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] != 0.0f)
 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_nnb  *nnb;
 655     int     i, j, k, j1, maxval;
 656
 657     /* Put all neighbors nearer than rmsdcut in the list */
 658     fprintf(stderr, "Making list of neighbors within cutoff ");
 659     snew(nnb, n1);
 660     for (i = 0; (i < n1); i++)
 661     {
 662         maxval = 0;
 663         k      = 0;
 664         /* put all neighbors within cut-off in list */
 665         for (j = 0; j < n1; j++)
 666         {
 667             if (mat[i][j] < rmsdcut)
 668             {
 669                 if (k >= maxval)
 670                 {
 671                     maxval += 10;
 672                     srenew(nnb[i].nb, maxval);
 673                 }
 674                 nnb[i].nb[k] = j;
 675                 k++;
 676             }
 677         }
 678         /* store nr of neighbors, we'll need that */
 679         nnb[i].nr = k;
 680         if (i%(1+n1/100) == 0)
 681         {
 682             fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
 683         }
 684     }
 685     fprintf(stderr, "%3d%%\n", 100);
 686
 687     /* sort neighbor list on number of neighbors, largest first */
 688     std::sort(nnb, nnb+n1, nrnb_comp);
 689
 690     if (debug)
 691     {
 692         dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
 693     }
 694
 695     /* turn first structure with all its neighbors (largest) into cluster
 696        remove them from pool of structures and repeat for all remaining */
 697     fprintf(stderr, "Finding clusters %4d", 0);
 698     /* cluster id's start at 1: */
 699     k = 1;
 700     while (nnb[0].nr)
 701     {
 702         /* set cluster id (k) for first item in neighborlist */
 703         for (j = 0; j < nnb[0].nr; j++)
 704         {
 705             clust->cl[nnb[0].nb[j]] = k;
 706         }
 707         /* mark as done */
 708         nnb[0].nr = 0;
 709         sfree(nnb[0].nb);
 710
 711         /* adjust number of neighbors for others, taking removals into account: */
 712         for (i = 1; i < n1 && nnb[i].nr; i++)
 713         {
 714             j1 = 0;
 715             for (j = 0; j < nnb[i].nr; j++)
 716             {
 717                 /* if this neighbor wasn't removed */
 718                 if (clust->cl[nnb[i].nb[j]] == 0)
 719                 {
 720                     /* shift the rest (j1<=j) */
 721                     nnb[i].nb[j1] = nnb[i].nb[j];
 722                     /* next */
 723                     j1++;
 724                 }
 725             }
 726             /* now j1 is the new number of neighbors */
 727             nnb[i].nr = j1;
 728         }
 729         /* sort again on nnb[].nr, because we have new # neighbors: */
 730         /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
 731         std::sort(nnb, nnb+i, nrnb_comp);
 732
 733         fprintf(stderr, "\b\b\b\b%4d", k);
 734         /* new cluster id */
 735         k++;
 736     }
 737     fprintf(stderr, "\n");
 738     sfree(nnb);
 739     if (debug)
 740     {
 741         fprintf(debug, "Clusters (%d):\n", k);
 742         for (i = 0; i < n1; i++)
 743         {
 744             fprintf(debug, " %3d", clust->cl[i]);
 745         }
 746         fprintf(debug, "\n");
 747     }
 748
 749     clust->ncl = k-1;
 750 }
 751
 752 static rvec **read_whole_trj(const char *fn, int isize, const int index[], int skip,
 753                              int *nframe, real **time,  matrix  **boxes, int **frameindices, const gmx_output_env_t *oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
 754 {
 755     rvec       **xx, *x;
 756     matrix       box;
 757     real         t;
 758     int          i, j, max_nf;
 759     int          natom;
 760     t_trxstatus *status;
 761
 762
 763     max_nf = 0;
 764     xx     = nullptr;
 765     *time  = nullptr;
 766     natom  = read_first_x(oenv, &status, fn, &t, &x, box);
 767     i      = 0;
 768     int clusterIndex = 0;
 769     do
 770     {
 771         if (bPBC)
 772         {
 773             gmx_rmpbc(gpbc, natom, box, x);
 774         }
 775         if (clusterIndex >= max_nf)
 776         {
 777             max_nf += 10;
 778             srenew(xx, max_nf);
 779             srenew(*time, max_nf);
 780             srenew(*boxes, max_nf);
 781             srenew(*frameindices, max_nf);
 782         }
 783         if ((i % skip) == 0)
 784         {
 785             snew(xx[clusterIndex], isize);
 786             /* Store only the interesting atoms */
 787             for (j = 0; (j < isize); j++)
 788             {
 789                 copy_rvec(x[index[j]], xx[clusterIndex][j]);
 790             }
 791             (*time)[clusterIndex] = t;
 792             copy_mat(box, (*boxes)[clusterIndex]);
 793             (*frameindices)[clusterIndex] = nframes_read(status);
 794             clusterIndex++;
 795         }
 796         i++;
 797     }
 798     while (read_next_x(oenv, status, &t, x, box));
 799     fprintf(stderr, "Allocated %zu bytes for frames\n",
 800             (max_nf*isize*sizeof(**xx)));
 801     fprintf(stderr, "Read %d frames from trajectory %s\n", clusterIndex, fn);
 802     *nframe = clusterIndex;
 803     sfree(x);
 804
 805     return xx;
 806 }
 807
 808 static int plot_clusters(int nf, real **mat, t_clusters *clust,
 809                          int minstruct)
 810 {
 811     int  i, j, ncluster, ci;
 812     int *cl_id, *nstruct, *strind;
 813
 814     snew(cl_id, nf);
 815     snew(nstruct, nf);
 816     snew(strind, nf);
 817     for (i = 0; i < nf; i++)
 818     {
 819         strind[i] = 0;
 820         cl_id[i]  = clust->cl[i];
 821         nstruct[cl_id[i]]++;
 822     }
 823     ncluster = 0;
 824     for (i = 0; i < nf; i++)
 825     {
 826         if (nstruct[i] >= minstruct)
 827         {
 828             ncluster++;
 829             for (j = 0; (j < nf); j++)
 830             {
 831                 if (cl_id[j] == i)
 832                 {
 833                     strind[j] = ncluster;
 834                 }
 835             }
 836         }
 837     }
 838     ncluster++;
 839     fprintf(stderr, "There are %d clusters with at least %d conformations\n",
 840             ncluster, minstruct);
 841
 842     for (i = 0; (i < nf); i++)
 843     {
 844         ci = cl_id[i];
 845         for (j = 0; j < i; j++)
 846         {
 847             if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
 848             {
 849                 /* color different clusters with different colors, as long as
 850                    we don't run out of colors */
 851                 mat[i][j] = strind[i];
 852             }
 853             else
 854             {
 855                 mat[i][j] = 0;
 856             }
 857         }
 858     }
 859     sfree(strind);
 860     sfree(nstruct);
 861     sfree(cl_id);
 862
 863     return ncluster;
 864 }
 865
 866 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
 867 {
 868     int i, j;
 869
 870     for (i = 0; i < nf; i++)
 871     {
 872         for (j = 0; j < i; j++)
 873         {
 874             if (clust->cl[i] == clust->cl[j])
 875             {
 876                 mat[i][j] = val;
 877             }
 878             else
 879             {
 880                 mat[i][j] = 0;
 881             }
 882         }
 883     }
 884 }
 885
 886 static char *parse_filename(const char *fn, int maxnr)
 887 {
 888     int         i;
 889     char       *fnout;
 890     const char *ext;
 891     char        buf[STRLEN];
 892
 893     if (std::strchr(fn, '%'))
 894     {
 895         gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
 896     }
 897     /* number of digits needed in numbering */
 898     i = static_cast<int>((std::log(static_cast<real>(maxnr))/std::log(10.0)) + 1);
 899     /* split fn and ext */
 900     ext = std::strrchr(fn, '.');
 901     if (!ext)
 902     {
 903         gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
 904     }
 905     ext++;
 906     /* insert e.g. '%03d' between fn and ext */
 907     sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
 908     snew(fnout, std::strlen(buf)+1);
 909     std::strcpy(fnout, buf);
 910
 911     return fnout;
 912 }
 913
 914 static void ana_trans(t_clusters *clust, int nf,
 915                       const char *transfn, const char *ntransfn, FILE *log,
 916                       t_rgb rlo, t_rgb rhi, const gmx_output_env_t *oenv)
 917 {
 918     FILE  *fp;
 919     real **trans, *axis;
 920     int   *ntrans;
 921     int    i, ntranst, maxtrans;
 922     char   buf[STRLEN];
 923
 924     snew(ntrans, clust->ncl);
 925     snew(trans, clust->ncl);
 926     snew(axis, clust->ncl);
 927     for (i = 0; i < clust->ncl; i++)
 928     {
 929         axis[i] = i+1;
 930         snew(trans[i], clust->ncl);
 931     }
 932     ntranst  = 0;
 933     maxtrans = 0;
 934     for (i = 1; i < nf; i++)
 935     {
 936         if (clust->cl[i] != clust->cl[i-1])
 937         {
 938             ntranst++;
 939             ntrans[clust->cl[i-1]-1]++;
 940             ntrans[clust->cl[i]-1]++;
 941             trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
 942             maxtrans = static_cast<int>(std::max(static_cast<real>(maxtrans), trans[clust->cl[i]-1][clust->cl[i-1]-1]));
 943         }
 944     }
 945     ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
 946                 "max %d between two specific clusters\n", ntranst, maxtrans);
 947     if (transfn)
 948     {
 949         fp = gmx_ffopen(transfn, "w");
 950         i  = std::min(maxtrans+1, 80);
 951         write_xpm(fp, 0, "Cluster Transitions", "# transitions",
 952                   "from cluster", "to cluster",
 953                   clust->ncl, clust->ncl, axis, axis, trans,
 954                   0, maxtrans, rlo, rhi, &i);
 955         gmx_ffclose(fp);
 956     }
 957     if (ntransfn)
 958     {
 959         fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
 960                       oenv);
 961         for (i = 0; i < clust->ncl; i++)
 962         {
 963             fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
 964         }
 965         xvgrclose(fp);
 966     }
 967     sfree(ntrans);
 968     for (i = 0; i < clust->ncl; i++)
 969     {
 970         sfree(trans[i]);
 971     }
 972     sfree(trans);
 973     sfree(axis);
 974 }
 975
 976 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
 977                              int natom, t_atoms *atoms, rvec *xtps,
 978                              real *mass, rvec **xx, real *time,
 979                              matrix *boxes, int *frameindices,
 980                              int ifsize, int *fitidx,
 981                              int iosize, int *outidx,
 982                              const char *trxfn, const char *sizefn,
 983                              const char *transfn, const char *ntransfn,
 984                              const char *clustidfn, const char *clustndxfn, gmx_bool bAverage,
 985                              int write_ncl, int write_nst, real rmsmin,
 986                              gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
 987                              const gmx_output_env_t *oenv)
 988 {
 989     FILE        *size_fp = nullptr;
 990     FILE        *ndxfn   = nullptr;
 991     char         buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
 992     t_trxstatus *trxout  = nullptr;
 993     t_trxstatus *trxsout = nullptr;
 994     int          i, i1, cl, nstr, *structure, first = 0, midstr;
 995     gmx_bool    *bWrite = nullptr;
 996     real         r, clrmsd, midrmsd;
 997     rvec        *xav = nullptr;
 998     matrix       zerobox;
 999
1000     clear_mat(zerobox);
1001
1002     ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1003     trxsfn = nullptr;
1004     if (trxfn)
1005     {
1006         /* do we write all structures? */
1007         if (write_ncl)
1008         {
1009             trxsfn = parse_filename(trxfn, std::max(write_ncl, clust->ncl));
1010             snew(bWrite, nf);
1011         }
1012         ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1013                     bAverage ? "average" : "middle", trxfn);
1014         if (write_ncl)
1015         {
1016             /* find out what we want to tell the user:
1017                Writing [all structures|structures with rmsd > %g] for
1018                {all|first %d} clusters {with more than %d structures} to %s     */
1019             if (rmsmin > 0.0)
1020             {
1021                 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1022             }
1023             else
1024             {
1025                 sprintf(buf1, "all structures");
1026             }
1027             buf2[0] = buf3[0] = '\0';
1028             if (write_ncl >= clust->ncl)
1029             {
1030                 if (write_nst == 0)
1031                 {
1032                     sprintf(buf2, "all ");
1033                 }
1034             }
1035             else
1036             {
1037                 sprintf(buf2, "the first %d ", write_ncl);
1038             }
1039             if (write_nst)
1040             {
1041                 sprintf(buf3, " with more than %d structures", write_nst);
1042             }
1043             sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1044             ffprintf(stderr, log, buf);
1045         }
1046
1047         /* Prepare a reference structure for the orientation of the clusters  */
1048         if (bFit)
1049         {
1050             reset_x(ifsize, fitidx, natom, nullptr, xtps, mass);
1051         }
1052         trxout = open_trx(trxfn, "w");
1053         /* Calculate the average structure in each cluster,               *
1054          * all structures are fitted to the first struture of the cluster */
1055         snew(xav, natom);
1056         GMX_ASSERT(xav, "");
1057     }
1058
1059     if (transfn || ntransfn)
1060     {
1061         ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1062     }
1063
1064     if (clustidfn)
1065     {
1066         FILE *fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1067         if (output_env_get_print_xvgr_codes(oenv))
1068         {
1069             fprintf(fp, "@    s0 symbol 2\n");
1070             fprintf(fp, "@    s0 symbol size 0.2\n");
1071             fprintf(fp, "@    s0 linestyle 0\n");
1072         }
1073         for (i = 0; i < nf; i++)
1074         {
1075             fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1076         }
1077         xvgrclose(fp);
1078     }
1079     if (sizefn)
1080     {
1081         size_fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1082         if (output_env_get_print_xvgr_codes(oenv))
1083         {
1084             fprintf(size_fp, "@g%d type %s\n", 0, "bar");
1085         }
1086     }
1087     if (clustndxfn && frameindices)
1088     {
1089         ndxfn = gmx_ffopen(clustndxfn, "w");
1090     }
1091
1092     snew(structure, nf);
1093     fprintf(log, "\n%3s | %3s  %4s | %6s %4s | cluster members\n",
1094             "cl.", "#st", "rmsd", "middle", "rmsd");
1095     for (cl = 1; cl <= clust->ncl; cl++)
1096     {
1097         /* prepare structures (fit, middle, average) */
1098         if (xav)
1099         {
1100             for (i = 0; i < natom; i++)
1101             {
1102                 clear_rvec(xav[i]);
1103             }
1104         }
1105         nstr = 0;
1106         for (i1 = 0; i1 < nf; i1++)
1107         {
1108             if (clust->cl[i1] == cl)
1109             {
1110                 structure[nstr] = i1;
1111                 nstr++;
1112                 if (trxfn && (bAverage || write_ncl) )
1113                 {
1114                     if (bFit)
1115                     {
1116                         reset_x(ifsize, fitidx, natom, nullptr, xx[i1], mass);
1117                     }
1118                     if (nstr == 1)
1119                     {
1120                         first = i1;
1121                     }
1122                     else if (bFit)
1123                     {
1124                         do_fit(natom, mass, xx[first], xx[i1]);
1125                     }
1126                     if (xav)
1127                     {
1128                         for (i = 0; i < natom; i++)
1129                         {
1130                             rvec_inc(xav[i], xx[i1][i]);
1131                         }
1132                     }
1133                 }
1134             }
1135         }
1136         if (sizefn)
1137         {
1138             fprintf(size_fp, "%8d %8d\n", cl, nstr);
1139         }
1140         if (ndxfn)
1141         {
1142             fprintf(ndxfn, "[Cluster_%04d]\n", cl);
1143         }
1144         clrmsd  = 0;
1145         midstr  = 0;
1146         midrmsd = 10000;
1147         for (i1 = 0; i1 < nstr; i1++)
1148         {
1149             r = 0;
1150             if (nstr > 1)
1151             {
1152                 for (i = 0; i < nstr; i++)
1153                 {
1154                     if (i < i1)
1155                     {
1156                         r += rmsd[structure[i]][structure[i1]];
1157                     }
1158                     else
1159                     {
1160                         r += rmsd[structure[i1]][structure[i]];
1161                     }
1162                 }
1163                 r /= (nstr - 1);
1164             }
1165             if (r < midrmsd)
1166             {
1167                 midstr  = structure[i1];
1168                 midrmsd = r;
1169             }
1170             clrmsd += r;
1171         }
1172         clrmsd /= nstr;
1173
1174         /* dump cluster info to logfile */
1175         if (nstr > 1)
1176         {
1177             sprintf(buf1, "%6.3f", clrmsd);
1178             if (buf1[0] == '0')
1179             {
1180                 buf1[0] = ' ';
1181             }
1182             sprintf(buf2, "%5.3f", midrmsd);
1183             if (buf2[0] == '0')
1184             {
1185                 buf2[0] = ' ';
1186             }
1187         }
1188         else
1189         {
1190             sprintf(buf1, "%5s", "");
1191             sprintf(buf2, "%5s", "");
1192         }
1193         fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1194         for (i = 0; i < nstr; i++)
1195         {
1196             if ((i % 7 == 0) && i)
1197             {
1198                 sprintf(buf, "\n%3s | %3s  %4s | %6s %4s |", "", "", "", "", "");
1199                 if (ndxfn)
1200                 {
1201                     fprintf(ndxfn, "\n");
1202                 }
1203             }
1204             else
1205             {
1206                 buf[0] = '\0';
1207             }
1208             i1 = structure[i];
1209             fprintf(log, "%s %6g", buf, time[i1]);
1210             if (ndxfn)
1211             {
1212                 fprintf(ndxfn, " %6d", frameindices[i1] + 1);
1213             }
1214         }
1215         fprintf(log, "\n");
1216         if (ndxfn)
1217         {
1218             fprintf(ndxfn, "\n");
1219         }
1220
1221         /* write structures to trajectory file(s) */
1222         if (trxfn)
1223         {
1224             if (write_ncl)
1225             {
1226                 for (i = 0; i < nstr; i++)
1227                 {
1228                     bWrite[i] = FALSE;
1229                 }
1230             }
1231             if (cl < write_ncl+1 && nstr > write_nst)
1232             {
1233                 /* Dump all structures for this cluster */
1234                 /* generate numbered filename (there is a %d in trxfn!) */
1235                 sprintf(buf, trxsfn, cl);
1236                 trxsout = open_trx(buf, "w");
1237                 for (i = 0; i < nstr; i++)
1238                 {
1239                     bWrite[i] = TRUE;
1240                     if (rmsmin > 0.0)
1241                     {
1242                         for (i1 = 0; i1 < i && bWrite[i]; i1++)
1243                         {
1244                             if (bWrite[i1])
1245                             {
1246                                 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1247                             }
1248                         }
1249                     }
1250                     if (bWrite[i])
1251                     {
1252                         write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], boxes[structure[i]],
1253                                   xx[structure[i]], nullptr, nullptr);
1254                     }
1255                 }
1256                 close_trx(trxsout);
1257             }
1258             /* Dump the average structure for this cluster */
1259             if (bAverage)
1260             {
1261                 for (i = 0; i < natom; i++)
1262                 {
1263                     svmul(1.0/nstr, xav[i], xav[i]);
1264                 }
1265             }
1266             else
1267             {
1268                 for (i = 0; i < natom; i++)
1269                 {
1270                     copy_rvec(xx[midstr][i], xav[i]);
1271                 }
1272                 if (bFit)
1273                 {
1274                     reset_x(ifsize, fitidx, natom, nullptr, xav, mass);
1275                 }
1276             }
1277             if (bFit)
1278             {
1279                 do_fit(natom, mass, xtps, xav);
1280             }
1281             write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], boxes[midstr], xav, nullptr, nullptr);
1282         }
1283     }
1284     /* clean up */
1285     if (trxfn)
1286     {
1287         close_trx(trxout);
1288         sfree(xav);
1289         if (write_ncl)
1290         {
1291             sfree(bWrite);
1292         }
1293     }
1294     sfree(structure);
1295     if (trxsfn)
1296     {
1297         sfree(trxsfn);
1298     }
1299
1300     if (size_fp)
1301     {
1302         xvgrclose(size_fp);
1303     }
1304     if (ndxfn)
1305     {
1306         gmx_ffclose(ndxfn);
1307     }
1308 }
1309
1310 static void convert_mat(t_matrix *mat, t_mat *rms)
1311 {
1312     int i, j;
1313
1314     rms->n1 = mat->nx;
1315     matrix2real(mat, rms->mat);
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  = std::max(rms->maxrms, rms->mat[i][j]);
1323             if (j != i)
1324             {
1325                 rms->minrms = std::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         "[THISMODULE] can cluster structures using several different methods.",
1336         "Distances between structures can be determined from a trajectory",
1337         "or read from an [REF].xpm[ref] 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         "[REF].xpm[ref] 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:",
1379         "",
1380         " * [TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1381         "   and a graphical depiction of the clusters in the lower right half",
1382         "   When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1383         "   when two structures are in the same cluster.",
1384         "   When [TT]-minstruct[tt] > 1 different colors will be used for each",
1385         "   cluster.",
1386         " * [TT]-g[tt] writes information on the options used and a detailed list",
1387         "   of all clusters and their members.",
1388         "",
1389
1390         "Additionally, a number of optional output files can be written:",
1391         "",
1392         " * [TT]-dist[tt] writes the RMSD distribution.",
1393         " * [TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1394         "   diagonalization.",
1395         " * [TT]-sz[tt] writes the cluster sizes.",
1396         " * [TT]-tr[tt] writes a matrix of the number transitions between",
1397         "   cluster pairs.",
1398         " * [TT]-ntr[tt] writes the total number of transitions to or from",
1399         "   each cluster.",
1400         " * [TT]-clid[tt] writes the cluster number as a function of time.",
1401         " * [TT]-clndx[tt] writes the frame numbers corresponding to the clusters to the",
1402         "   specified index file to be read into trjconv.",
1403         " * [TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1404         "   structure of each cluster or writes numbered files with cluster members",
1405         "   for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1406         "   [TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1407         "   structure with the smallest average RMSD from all other structures",
1408         "   of the cluster.",
1409     };
1410
1411     FILE              *fp, *log;
1412     int                nf   = 0, i, i1, i2, j;
1413     int64_t            nrms = 0;
1414
1415     matrix             box;
1416     matrix            *boxes = nullptr;
1417     rvec              *xtps, *usextps, *x1, **xx = nullptr;
1418     const char        *fn, *trx_out_fn;
1419     t_clusters         clust;
1420     t_mat             *rms, *orig = nullptr;
1421     real              *eigenvalues;
1422     t_topology         top;
1423     int                ePBC;
1424     t_atoms            useatoms;
1425     real              *eigenvectors;
1426
1427     int                isize = 0, ifsize = 0, iosize = 0;
1428     int               *index = nullptr, *fitidx = nullptr, *outidx = nullptr, *frameindices = nullptr;
1429     char              *grpname;
1430     real               rmsd, **d1, **d2, *time = nullptr, time_invfac, *mass = nullptr;
1431     char               buf[STRLEN], buf1[80];
1432     gmx_bool           bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1433
1434     int                method, ncluster = 0;
1435     static const char *methodname[] = {
1436         nullptr, "linkage", "jarvis-patrick", "monte-carlo",
1437         "diagonalization", "gromos", nullptr
1438     };
1439     enum {
1440         m_null, m_linkage, m_jarvis_patrick,
1441         m_monte_carlo, m_diagonalize, m_gromos, m_nr
1442     };
1443     /* Set colors for plotting: white = zero RMS, black = maximum */
1444     static t_rgb      rlo_top  = { 1.0, 1.0, 1.0 };
1445     static t_rgb      rhi_top  = { 0.0, 0.0, 0.0 };
1446     static t_rgb      rlo_bot  = { 1.0, 1.0, 1.0 };
1447     static t_rgb      rhi_bot  = { 0.0, 0.0, 1.0 };
1448     static int        nlevels  = 40, skip = 1;
1449     static real       scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1450     gmx_bool          bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1451     static int        niter    = 10000, nrandom = 0, seed = 0, write_ncl = 0, write_nst = 1, minstruct = 1;
1452     static real       kT       = 1e-3;
1453     static int        M        = 10, P = 3;
1454     gmx_output_env_t *oenv;
1455     gmx_rmpbc_t       gpbc = nullptr;
1456
1457     t_pargs           pa[] = {
1458         { "-dista", FALSE, etBOOL, {&bRMSdist},
1459           "Use RMSD of distances instead of RMS deviation" },
1460         { "-nlevels", FALSE, etINT,  {&nlevels},
1461           "Discretize RMSD matrix in this number of levels" },
1462         { "-cutoff", FALSE, etREAL, {&rmsdcut},
1463           "RMSD cut-off (nm) for two structures to be neighbor" },
1464         { "-fit",   FALSE, etBOOL, {&bFit},
1465           "Use least squares fitting before RMSD calculation" },
1466         { "-max",   FALSE, etREAL, {&scalemax},
1467           "Maximum level in RMSD matrix" },
1468         { "-skip",  FALSE, etINT,  {&skip},
1469           "Only analyze every nr-th frame" },
1470         { "-av",    FALSE, etBOOL, {&bAverage},
1471           "Write average instead of middle structure for each cluster" },
1472         { "-wcl",   FALSE, etINT,  {&write_ncl},
1473           "Write the structures for this number of clusters to numbered files" },
1474         { "-nst",   FALSE, etINT,  {&write_nst},
1475           "Only write all structures if more than this number of structures per cluster" },
1476         { "-rmsmin", FALSE, etREAL, {&rmsmin},
1477           "minimum rms difference with rest of cluster for writing structures" },
1478         { "-method", FALSE, etENUM, {methodname},
1479           "Method for cluster determination" },
1480         { "-minstruct", FALSE, etINT, {&minstruct},
1481           "Minimum number of structures in cluster for coloring in the [REF].xpm[ref] file" },
1482         { "-binary", FALSE, etBOOL, {&bBinary},
1483           "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1484           "is given by [TT]-cutoff[tt]" },
1485         { "-M",     FALSE, etINT,  {&M},
1486           "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1487           "0 is use cutoff" },
1488         { "-P",     FALSE, etINT,  {&P},
1489           "Number of identical nearest neighbors required to form a cluster" },
1490         { "-seed",  FALSE, etINT,  {&seed},
1491           "Random number seed for Monte Carlo clustering algorithm (0 means generate)" },
1492         { "-niter", FALSE, etINT,  {&niter},
1493           "Number of iterations for MC" },
1494         { "-nrandom", FALSE, etINT,  {&nrandom},
1495           "The first iterations for MC may be done complete random, to shuffle the frames" },
1496         { "-kT",    FALSE, etREAL, {&kT},
1497           "Boltzmann weighting factor for Monte Carlo optimization "
1498           "(zero turns off uphill steps)" },
1499         { "-pbc", FALSE, etBOOL,
1500           { &bPBC }, "PBC check" }
1501     };
1502     t_filenm          fnm[] = {
1503         { efTRX, "-f",     nullptr,        ffOPTRD },
1504         { efTPS, "-s",     nullptr,        ffREAD },
1505         { efNDX, nullptr,     nullptr,        ffOPTRD },
1506         { efXPM, "-dm",   "rmsd",       ffOPTRD },
1507         { efXPM, "-om",   "rmsd-raw",   ffWRITE },
1508         { efXPM, "-o",    "rmsd-clust", ffWRITE },
1509         { efLOG, "-g",    "cluster",    ffWRITE },
1510         { efXVG, "-dist", "rmsd-dist",  ffOPTWR },
1511         { efXVG, "-ev",   "rmsd-eig",   ffOPTWR },
1512         { efXVG, "-conv", "mc-conv",    ffOPTWR },
1513         { efXVG, "-sz",   "clust-size", ffOPTWR},
1514         { efXPM, "-tr",   "clust-trans", ffOPTWR},
1515         { efXVG, "-ntr",  "clust-trans", ffOPTWR},
1516         { efXVG, "-clid", "clust-id",   ffOPTWR},
1517         { efTRX, "-cl",   "clusters.pdb", ffOPTWR },
1518         { efNDX, "-clndx", "clusters.ndx", ffOPTWR }
1519     };
1520 #define NFILE asize(fnm)
1521
1522     if (!parse_common_args(&argc, argv,
1523                            PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
1524                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr,
1525                            &oenv))
1526     {
1527         return 0;
1528     }
1529
1530     /* parse options */
1531     bReadMat   = opt2bSet("-dm", NFILE, fnm);
1532     bReadTraj  = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1533     if (opt2parg_bSet("-av", asize(pa), pa) ||
1534         opt2parg_bSet("-wcl", asize(pa), pa) ||
1535         opt2parg_bSet("-nst", asize(pa), pa) ||
1536         opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1537         opt2bSet("-cl", NFILE, fnm) )
1538     {
1539         trx_out_fn = opt2fn("-cl", NFILE, fnm);
1540     }
1541     else
1542     {
1543         trx_out_fn = nullptr;
1544     }
1545     if (bReadMat && output_env_get_time_factor(oenv) != 1)
1546     {
1547         fprintf(stderr,
1548                 "\nWarning: assuming the time unit in %s is %s\n",
1549                 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv).c_str());
1550     }
1551     if (trx_out_fn && !bReadTraj)
1552     {
1553         fprintf(stderr, "\nWarning: "
1554                 "cannot write cluster structures without reading trajectory\n"
1555                 "         ignoring option -cl %s\n", trx_out_fn);
1556     }
1557
1558     method = 1;
1559     while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1560     {
1561         method++;
1562     }
1563     if (method == m_nr)
1564     {
1565         gmx_fatal(FARGS, "Invalid method");
1566     }
1567
1568     bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1569                 method == m_gromos );
1570
1571     /* Open log file */
1572     log = ftp2FILE(efLOG, NFILE, fnm, "w");
1573
1574     fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1575     fprintf(log, "Using %s method for clustering\n", methodname[0]);
1576
1577     /* check input and write parameters to log file */
1578     bUseRmsdCut = FALSE;
1579     if (method == m_jarvis_patrick)
1580     {
1581         bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1582         if ((M < 0) || (M == 1))
1583         {
1584             gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1585         }
1586         if (M < 2)
1587         {
1588             sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1589             bUseRmsdCut = TRUE;
1590         }
1591         else
1592         {
1593             if (P >= M)
1594             {
1595                 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1596             }
1597             if (bJP_RMSD)
1598             {
1599                 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1600                 bUseRmsdCut = TRUE;
1601             }
1602             else
1603             {
1604                 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1605             }
1606         }
1607         ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1608     }
1609     else /* method != m_jarvis */
1610     {
1611         bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1612     }
1613     if (bUseRmsdCut && method != m_jarvis_patrick)
1614     {
1615         fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1616     }
1617     if (method == m_monte_carlo)
1618     {
1619         fprintf(log, "Using %d iterations\n", niter);
1620     }
1621
1622     if (skip < 1)
1623     {
1624         gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1625     }
1626
1627     /* get input */
1628     if (bReadTraj)
1629     {
1630         /* don't read mass-database as masses (and top) are not used */
1631         read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtps, nullptr, box,
1632                       TRUE);
1633         if (bPBC)
1634         {
1635             gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1636         }
1637
1638         fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1639                 bReadMat ? "" : " and RMSD calculation");
1640         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1641                   1, &ifsize, &fitidx, &grpname);
1642         if (trx_out_fn)
1643         {
1644             fprintf(stderr, "\nSelect group for output:\n");
1645             get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1646                       1, &iosize, &outidx, &grpname);
1647             /* merge and convert both index groups: */
1648             /* first copy outidx to index. let outidx refer to elements in index */
1649             snew(index, iosize);
1650             isize = iosize;
1651             for (i = 0; i < iosize; i++)
1652             {
1653                 index[i]  = outidx[i];
1654                 outidx[i] = i;
1655             }
1656             /* now lookup elements from fitidx in index, add them if necessary
1657                and also let fitidx refer to elements in index */
1658             for (i = 0; i < ifsize; i++)
1659             {
1660                 j = 0;
1661                 while (j < isize && index[j] != fitidx[i])
1662                 {
1663                     j++;
1664                 }
1665                 if (j >= isize)
1666                 {
1667                     /* slow this way, but doesn't matter much */
1668                     isize++;
1669                     srenew(index, isize);
1670                 }
1671                 index[j]  = fitidx[i];
1672                 fitidx[i] = j;
1673             }
1674         }
1675         else /* !trx_out_fn */
1676         {
1677             isize = ifsize;
1678             snew(index, isize);
1679             for (i = 0; i < ifsize; i++)
1680             {
1681                 index[i]  = fitidx[i];
1682                 fitidx[i] = i;
1683             }
1684         }
1685     }
1686
1687     if (bReadTraj)
1688     {
1689         /* Loop over first coordinate file */
1690         fn = opt2fn("-f", NFILE, fnm);
1691
1692         xx = read_whole_trj(fn, isize, index, skip, &nf, &time, &boxes, &frameindices, oenv, bPBC, gpbc);
1693         output_env_conv_times(oenv, nf, time);
1694         if (!bRMSdist || bAnalyze)
1695         {
1696             /* Center all frames on zero */
1697             snew(mass, isize);
1698             for (i = 0; i < ifsize; i++)
1699             {
1700                 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1701             }
1702             if (bFit)
1703             {
1704                 for (i = 0; i < nf; i++)
1705                 {
1706                     reset_x(ifsize, fitidx, isize, nullptr, xx[i], mass);
1707                 }
1708             }
1709         }
1710         if (bPBC)
1711         {
1712             gmx_rmpbc_done(gpbc);
1713         }
1714     }
1715
1716     std::vector<t_matrix> readmat;
1717     if (bReadMat)
1718     {
1719         fprintf(stderr, "Reading rms distance matrix ");
1720         readmat = read_xpm_matrix(opt2fn("-dm", NFILE, fnm));
1721         fprintf(stderr, "\n");
1722         if (readmat[0].nx != readmat[0].ny)
1723         {
1724             gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1725                       readmat[0].nx, readmat[0].ny);
1726         }
1727         if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1728         {
1729             gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1730                       "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1731         }
1732
1733         nf = readmat[0].nx;
1734         sfree(time);
1735         time        = readmat[0].axis_x.data();
1736         time_invfac = output_env_get_time_invfactor(oenv);
1737         for (i = 0; i < nf; i++)
1738         {
1739             time[i] *= time_invfac;
1740         }
1741
1742         rms = init_mat(readmat[0].nx, method == m_diagonalize);
1743         convert_mat(&(readmat[0]), rms);
1744
1745         nlevels = gmx::ssize(readmat[0].map);
1746     }
1747     else   /* !bReadMat */
1748     {
1749         rms  = init_mat(nf, method == m_diagonalize);
1750         nrms = (static_cast<int64_t>(nf)*static_cast<int64_t>(nf-1))/2;
1751         if (!bRMSdist)
1752         {
1753             fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1754             /* Initialize work array */
1755             snew(x1, isize);
1756             for (i1 = 0; i1 < nf; i1++)
1757             {
1758                 for (i2 = i1+1; i2 < nf; i2++)
1759                 {
1760                     for (i = 0; i < isize; i++)
1761                     {
1762                         copy_rvec(xx[i1][i], x1[i]);
1763                     }
1764                     if (bFit)
1765                     {
1766                         do_fit(isize, mass, xx[i2], x1);
1767                     }
1768                     rmsd = rmsdev(isize, mass, xx[i2], x1);
1769                     set_mat_entry(rms, i1, i2, rmsd);
1770                 }
1771                 nrms -= nf-i1-1;
1772                 fprintf(stderr, "\r# RMSD calculations left: " "%" PRId64 "   ", nrms);
1773                 fflush(stderr);
1774             }
1775             sfree(x1);
1776         }
1777         else /* bRMSdist */
1778         {
1779             fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1780
1781             /* Initiate work arrays */
1782             snew(d1, isize);
1783             snew(d2, isize);
1784             for (i = 0; (i < isize); i++)
1785             {
1786                 snew(d1[i], isize);
1787                 snew(d2[i], isize);
1788             }
1789             for (i1 = 0; i1 < nf; i1++)
1790             {
1791                 calc_dist(isize, xx[i1], d1);
1792                 for (i2 = i1+1; (i2 < nf); i2++)
1793                 {
1794                     calc_dist(isize, xx[i2], d2);
1795                     set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1796                 }
1797                 nrms -= nf-i1-1;
1798                 fprintf(stderr, "\r# RMSD calculations left: " "%" PRId64 "   ", nrms);
1799                 fflush(stderr);
1800             }
1801             /* Clean up work arrays */
1802             for (i = 0; (i < isize); i++)
1803             {
1804                 sfree(d1[i]);
1805                 sfree(d2[i]);
1806             }
1807             sfree(d1);
1808             sfree(d2);
1809         }
1810         fprintf(stderr, "\n\n");
1811     }
1812     ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1813                 rms->minrms, rms->maxrms);
1814     ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1815     ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1816     ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1817     if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1818     {
1819         fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1820                 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1821     }
1822     if (bAnalyze && (rmsmin < rms->minrms) )
1823     {
1824         fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1825                 rmsmin, rms->minrms);
1826     }
1827     if (bAnalyze && (rmsmin > rmsdcut) )
1828     {
1829         fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1830                 rmsmin, rmsdcut);
1831     }
1832
1833     /* Plot the rmsd distribution */
1834     rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1835
1836     if (bBinary)
1837     {
1838         for (i1 = 0; (i1 < nf); i1++)
1839         {
1840             for (i2 = 0; (i2 < nf); i2++)
1841             {
1842                 if (rms->mat[i1][i2] < rmsdcut)
1843                 {
1844                     rms->mat[i1][i2] = 0;
1845                 }
1846                 else
1847                 {
1848                     rms->mat[i1][i2] = 1;
1849                 }
1850             }
1851         }
1852     }
1853
1854     snew(clust.cl, nf);
1855     switch (method)
1856     {
1857         case m_linkage:
1858             /* Now sort the matrix and write it out again */
1859             gather(rms, rmsdcut, &clust);
1860             break;
1861         case m_diagonalize:
1862             /* Do a diagonalization */
1863             snew(eigenvalues, nf);
1864             snew(eigenvectors, nf*nf);
1865             std::memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1866             eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1867             sfree(eigenvectors);
1868
1869             fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1870                           "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1871             for (i = 0; (i < nf); i++)
1872             {
1873                 fprintf(fp, "%10d  %10g\n", i, eigenvalues[i]);
1874             }
1875             xvgrclose(fp);
1876             break;
1877         case m_monte_carlo:
1878             orig     = init_mat(rms->nn, FALSE);
1879             orig->nn = rms->nn;
1880             copy_t_mat(orig, rms);
1881             mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1882                         opt2fn_null("-conv", NFILE, fnm), oenv);
1883             break;
1884         case m_jarvis_patrick:
1885             jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1886             break;
1887         case m_gromos:
1888             gromos(rms->nn, rms->mat, rmsdcut, &clust);
1889             break;
1890         default:
1891             gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1892     }
1893
1894     if (method == m_monte_carlo || method == m_diagonalize)
1895     {
1896         fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1897                 mat_energy(rms));
1898     }
1899
1900     if (bAnalyze)
1901     {
1902         if (minstruct > 1)
1903         {
1904             ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1905         }
1906         else
1907         {
1908             mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1909         }
1910         init_t_atoms(&useatoms, isize, FALSE);
1911         snew(usextps, isize);
1912         useatoms.resinfo = top.atoms.resinfo;
1913         for (i = 0; i < isize; i++)
1914         {
1915             useatoms.atomname[i]    = top.atoms.atomname[index[i]];
1916             useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1917             useatoms.nres           = std::max(useatoms.nres, useatoms.atom[i].resind+1);
1918             copy_rvec(xtps[index[i]], usextps[i]);
1919         }
1920         useatoms.nr = isize;
1921         analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time, boxes, frameindices,
1922                          ifsize, fitidx, iosize, outidx,
1923                          bReadTraj ? trx_out_fn : nullptr,
1924                          opt2fn_null("-sz", NFILE, fnm),
1925                          opt2fn_null("-tr", NFILE, fnm),
1926                          opt2fn_null("-ntr", NFILE, fnm),
1927                          opt2fn_null("-clid", NFILE, fnm),
1928                          opt2fn_null("-clndx", NFILE, fnm),
1929                          bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1930                          rlo_bot, rhi_bot, oenv);
1931         sfree(boxes);
1932         sfree(frameindices);
1933     }
1934     gmx_ffclose(log);
1935
1936     if (bBinary && !bAnalyze)
1937     {
1938         /* Make the clustering visible */
1939         for (i2 = 0; (i2 < nf); i2++)
1940         {
1941             for (i1 = i2+1; (i1 < nf); i1++)
1942             {
1943                 if (rms->mat[i1][i2] != 0.0f)
1944                 {
1945                     rms->mat[i1][i2] = rms->maxrms;
1946                 }
1947             }
1948         }
1949     }
1950
1951     fp = opt2FILE("-o", NFILE, fnm, "w");
1952     fprintf(stderr, "Writing rms distance/clustering matrix ");
1953     if (bReadMat)
1954     {
1955         write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1956                   readmat[0].label_y, nf, nf, readmat[0].axis_x.data(), readmat[0].axis_y.data(),
1957                   rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1958     }
1959     else
1960     {
1961         auto timeLabel = output_env_get_time_label(oenv);
1962         auto title     = gmx::formatString("RMS%sDeviation / Cluster Index",
1963                                            bRMSdist ? " Distance " : " ");
1964         if (minstruct > 1)
1965         {
1966             write_xpm_split(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1967                             nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1968                             rlo_top, rhi_top, 0.0, ncluster,
1969                             &ncluster, TRUE, rlo_bot, rhi_bot);
1970         }
1971         else
1972         {
1973             write_xpm(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1974                       nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1975                       rlo_top, rhi_top, &nlevels);
1976         }
1977     }
1978     fprintf(stderr, "\n");
1979     gmx_ffclose(fp);
1980     if (nullptr != orig)
1981     {
1982         fp = opt2FILE("-om", NFILE, fnm, "w");
1983         auto timeLabel = output_env_get_time_label(oenv);
1984         auto title     = gmx::formatString("RMS%sDeviation", bRMSdist ? " Distance " : " ");
1985         write_xpm(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1986                   nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1987                   rlo_top, rhi_top, &nlevels);
1988         gmx_ffclose(fp);
1989         done_mat(&orig);
1990         sfree(orig);
1991     }
1992     /* now show what we've done */
1993     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1994     do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1995     if (method == m_diagonalize)
1996     {
1997         do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1998     }
1999     do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
2000     if (bAnalyze)
2001     {
2002         do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
2003         do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
2004         do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
2005     }
2006     do_view(oenv, opt2fn_null("-conv", NFILE, fnm), nullptr);
2007
2008     return 0;
2009 }