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