2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
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"
70 /* print to two file pointers at once (i.e. stderr and log) */
72 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
74 fprintf(fp1, "%s", buf);
75 fprintf(fp2, "%s", buf);
78 /* just print a prepared buffer to fp1 and fp2 */
80 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
82 lo_ffprintf(fp1, fp2, buf);
85 /* prepare buffer with one argument, then print to fp1 and fp2 */
87 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
89 sprintf(buf, fmt, arg);
90 lo_ffprintf(fp1, fp2, buf);
93 /* prepare buffer with one argument, then print to fp1 and fp2 */
95 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
97 sprintf(buf, fmt, arg);
98 lo_ffprintf(fp1, fp2, buf);
101 /* prepare buffer with one argument, then print to fp1 and fp2 */
103 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
105 sprintf(buf, fmt, arg);
106 lo_ffprintf(fp1, fp2, buf);
109 /* prepare buffer with two arguments, then print to fp1 and fp2 */
111 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
113 sprintf(buf, fmt, arg1, arg2);
114 lo_ffprintf(fp1, fp2, buf);
117 /* prepare buffer with two arguments, then print to fp1 and fp2 */
119 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
121 sprintf(buf, fmt, arg1, arg2);
122 lo_ffprintf(fp1, fp2, buf);
125 /* prepare buffer with two arguments, then print to fp1 and fp2 */
127 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
129 sprintf(buf, fmt, arg1, arg2);
130 lo_ffprintf(fp1, fp2, buf);
143 void cp_index(int nn, int from[], int to[])
147 for (i = 0; (i < nn); i++)
153 void mc_optimize(FILE *log, t_mat *m, real *time,
154 int maxiter, int nrandom,
156 const char *conv, gmx_output_env_t *oenv)
159 real ecur, enext, emin, prob, enorm;
160 int i, j, iswap, jswap, nn, nuphill = 0;
165 seed = static_cast<int>(gmx::makeRandomSeed());
167 gmx::DefaultRandomEngine rng(seed);
171 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
174 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
175 printf("by reordering the frames to minimize the path between the two structures\n");
176 printf("that have the largest pairwise RMSD.\n");
177 printf("Using random seed %d.\n", seed);
180 enorm = m->mat[0][0];
181 for (i = 0; (i < m->n1); i++)
183 for (j = 0; (j < m->nn); j++)
185 if (m->mat[i][j] > enorm)
187 enorm = m->mat[i][j];
193 if ((iswap == -1) || (jswap == -1))
195 fprintf(stderr, "Matrix contains identical values in all fields\n");
198 swap_rows(m, 0, iswap);
199 swap_rows(m, m->n1-1, jswap);
200 emin = ecur = mat_energy(m);
201 printf("Largest distance %g between %d and %d. Energy: %g.\n",
202 enorm, iswap, jswap, emin);
206 /* Initiate and store global minimum */
207 minimum = init_mat(nn, m->b1D);
209 copy_t_mat(minimum, m);
213 fp = xvgropen(conv, "Convergence of the MC optimization",
214 "Energy", "Step", oenv);
217 gmx::UniformIntDistribution<int> intDistNN(1, nn-2); // [1,nn-2]
218 gmx::UniformRealDistribution<real> realDistOne; // [0,1)
220 for (i = 0; (i < maxiter); i++)
222 /* Generate new swapping candidates */
225 iswap = intDistNN(rng);
226 jswap = intDistNN(rng);
228 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
230 /* Apply swap and compute energy */
231 swap_rows(m, iswap, jswap);
232 enext = mat_energy(m);
234 /* Compute probability */
236 if ((enext < ecur) || (i < nrandom))
241 /* Store global minimum */
242 copy_t_mat(minimum, m);
248 /* Try Monte Carlo step */
249 prob = std::exp(-(enext-ecur)/(enorm*kT));
252 if (prob == 1 || realDistOne(rng) < prob)
259 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
260 i, iswap, jswap, enext, prob);
263 fprintf(fp, "%6d %10g\n", i, enext);
269 swap_rows(m, jswap, iswap);
272 fprintf(log, "%d uphill steps were taken during optimization\n",
275 /* Now swap the matrix to get it into global minimum mode */
276 copy_t_mat(m, minimum);
278 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
279 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
280 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
281 for (i = 0; (i < m->nn); i++)
283 fprintf(log, "%10g %5d %10g\n",
286 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
295 static void calc_dist(int nind, rvec x[], real **d)
301 for (i = 0; (i < nind-1); i++)
304 for (j = i+1; (j < nind); j++)
306 /* Should use pbc_dx when analysing multiple molecueles,
307 * but the box is not stored for every frame.
309 rvec_sub(xi, x[j], dx);
315 static real rms_dist(int isize, real **d, real **d_r)
321 for (i = 0; (i < isize-1); i++)
323 for (j = i+1; (j < isize); j++)
325 r = d[i][j]-d_r[i][j];
329 r2 /= (isize*(isize-1))/2;
331 return std::sqrt(r2);
334 static bool rms_dist_comp(const t_dist &a, const t_dist &b)
336 return a.dist < b.dist;
339 static bool clust_id_comp(const t_clustid &a, const t_clustid &b)
341 return a.clust < b.clust;
344 static bool nrnb_comp(const t_nnb &a, const t_nnb &b)
346 /* return b<a, we want highest first */
350 void gather(t_mat *m, real cutoff, t_clusters *clust)
354 int i, j, k, nn, cid, n1, diff;
357 /* First we sort the entries in the RMSD matrix */
361 for (i = k = 0; (i < n1); i++)
363 for (j = i+1; (j < n1); j++, k++)
367 d[k].dist = m->mat[i][j];
372 gmx_incons("gather algortihm");
374 std::sort(d, d+nn, rms_dist_comp);
376 /* Now we make a cluster index for all of the conformations */
379 /* Now we check the closest structures, and equalize their cluster numbers */
380 fprintf(stderr, "Linking structures ");
383 fprintf(stderr, "*");
385 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
387 diff = c[d[k].j].clust - c[d[k].i].clust;
393 c[d[k].j].clust = c[d[k].i].clust;
397 c[d[k].i].clust = c[d[k].j].clust;
403 fprintf(stderr, "\nSorting and renumbering clusters\n");
404 /* Sort on cluster number */
405 std::sort(c, c+n1, clust_id_comp);
407 /* Renumber clusters */
409 for (k = 1; k < n1; k++)
411 if (c[k].clust != c[k-1].clust)
424 for (k = 0; (k < n1); k++)
426 fprintf(debug, "Cluster index for conformation %d: %d\n",
427 c[k].conf, c[k].clust);
431 for (k = 0; k < n1; k++)
433 clust->cl[c[k].conf] = c[k].clust;
440 gmx_bool jp_same(int **nnb, int i, int j, int P)
446 for (k = 0; nnb[i][k] >= 0; k++)
448 bIn = bIn || (nnb[i][k] == j);
456 for (k = 0; nnb[j][k] >= 0; k++)
458 bIn = bIn || (nnb[j][k] == i);
466 for (ii = 0; nnb[i][ii] >= 0; ii++)
468 for (jj = 0; nnb[j][jj] >= 0; jj++)
470 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
480 static void jarvis_patrick(int n1, real **mat, int M, int P,
481 real rmsdcut, t_clusters *clust)
486 int i, j, k, cid, diff, maxval;
488 real **mcpy = nullptr;
495 /* First we sort the entries in the RMSD matrix row by row.
496 * This gives us the nearest neighbor list.
500 for (i = 0; (i < n1); i++)
502 for (j = 0; (j < n1); j++)
505 row[j].dist = mat[i][j];
507 std::sort(row, row+n1, rms_dist_comp);
510 /* Put the M nearest neighbors in the list */
512 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
516 nnb[i][k] = row[j].j;
524 /* Put all neighbors nearer than rmsdcut in the list */
527 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
534 srenew(nnb[i], maxval);
536 nnb[i][k] = row[j].j;
542 srenew(nnb[i], maxval+1);
550 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
551 for (i = 0; (i < n1); i++)
553 fprintf(debug, "i:%5d nbs:", i);
554 for (j = 0; nnb[i][j] >= 0; j++)
556 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
558 fprintf(debug, "\n");
563 fprintf(stderr, "Linking structures ");
564 /* Use mcpy for temporary storage of booleans */
565 mcpy = mk_matrix(n1, n1, FALSE);
566 for (i = 0; i < n1; i++)
568 for (j = i+1; j < n1; j++)
570 mcpy[i][j] = jp_same(nnb, i, j, P);
575 fprintf(stderr, "*");
577 for (i = 0; i < n1; i++)
579 for (j = i+1; j < n1; j++)
583 diff = c[j].clust - c[i].clust;
589 c[j].clust = c[i].clust;
593 c[i].clust = c[j].clust;
602 fprintf(stderr, "\nSorting and renumbering clusters\n");
603 /* Sort on cluster number */
604 std::sort(c, c+n1, clust_id_comp);
606 /* Renumber clusters */
608 for (k = 1; k < n1; k++)
610 if (c[k].clust != c[k-1].clust)
622 for (k = 0; k < n1; k++)
624 clust->cl[c[k].conf] = c[k].clust;
628 for (k = 0; (k < n1); k++)
630 fprintf(debug, "Cluster index for conformation %d: %d\n",
631 c[k].conf, c[k].clust);
635 done_matrix(n1, &mcpy);
638 for (i = 0; (i < n1); i++)
645 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
649 /* dump neighbor list */
650 fprintf(fp, "%s", title);
651 for (i = 0; (i < n1); i++)
653 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
654 for (j = 0; j < nnb[i].nr; j++)
656 fprintf(fp, "%5d", nnb[i].nb[j]);
662 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
666 int i, j, k, j1, maxval;
668 /* Put all neighbors nearer than rmsdcut in the list */
669 fprintf(stderr, "Making list of neighbors within cutoff ");
672 for (i = 0; (i < n1); i++)
676 /* put all neighbors within cut-off in list */
677 for (j = 0; j < n1; j++)
679 if (mat[i][j] < rmsdcut)
684 srenew(nnb[i].nb, maxval);
690 /* store nr of neighbors, we'll need that */
692 if (i%(1+n1/100) == 0)
694 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
697 fprintf(stderr, "%3d%%\n", 100);
700 /* sort neighbor list on number of neighbors, largest first */
701 std::sort(nnb, nnb+n1, nrnb_comp);
705 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
708 /* turn first structure with all its neighbors (largest) into cluster
709 remove them from pool of structures and repeat for all remaining */
710 fprintf(stderr, "Finding clusters %4d", 0);
711 /* cluster id's start at 1: */
715 /* set cluster id (k) for first item in neighborlist */
716 for (j = 0; j < nnb[0].nr; j++)
718 clust->cl[nnb[0].nb[j]] = k;
724 /* adjust number of neighbors for others, taking removals into account: */
725 for (i = 1; i < n1 && nnb[i].nr; i++)
728 for (j = 0; j < nnb[i].nr; j++)
730 /* if this neighbor wasn't removed */
731 if (clust->cl[nnb[i].nb[j]] == 0)
733 /* shift the rest (j1<=j) */
734 nnb[i].nb[j1] = nnb[i].nb[j];
739 /* now j1 is the new number of neighbors */
742 /* sort again on nnb[].nr, because we have new # neighbors: */
743 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
744 std::sort(nnb, nnb+i, nrnb_comp);
746 fprintf(stderr, "\b\b\b\b%4d", k);
750 fprintf(stderr, "\n");
754 fprintf(debug, "Clusters (%d):\n", k);
755 for (i = 0; i < n1; i++)
757 fprintf(debug, " %3d", clust->cl[i]);
759 fprintf(debug, "\n");
765 rvec **read_whole_trj(const char *fn, int isize, int index[], int skip,
766 int *nframe, real **time, const gmx_output_env_t *oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
771 int i, i0, j, max_nf;
779 natom = read_first_x(oenv, &status, fn, &t, &x, box);
786 gmx_rmpbc(gpbc, natom, box, x);
792 srenew(*time, max_nf);
797 /* Store only the interesting atoms */
798 for (j = 0; (j < isize); j++)
800 copy_rvec(x[index[j]], xx[i0][j]);
807 while (read_next_x(oenv, status, &t, x, box));
808 fprintf(stderr, "Allocated %lu bytes for frames\n",
809 (unsigned long) (max_nf*isize*sizeof(**xx)));
810 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
817 static int plot_clusters(int nf, real **mat, t_clusters *clust,
820 int i, j, ncluster, ci;
821 int *cl_id, *nstruct, *strind;
826 for (i = 0; i < nf; i++)
829 cl_id[i] = clust->cl[i];
833 for (i = 0; i < nf; i++)
835 if (nstruct[i] >= minstruct)
838 for (j = 0; (j < nf); j++)
842 strind[j] = ncluster;
848 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
849 ncluster, minstruct);
851 for (i = 0; (i < nf); i++)
854 for (j = 0; j < i; j++)
856 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
858 /* color different clusters with different colors, as long as
859 we don't run out of colors */
860 mat[i][j] = strind[i];
875 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
879 for (i = 0; i < nf; i++)
881 for (j = 0; j < i; j++)
883 if (clust->cl[i] == clust->cl[j])
895 static char *parse_filename(const char *fn, int maxnr)
902 if (std::strchr(fn, '%'))
904 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
906 /* number of digits needed in numbering */
907 i = static_cast<int>((std::log(static_cast<real>(maxnr))/std::log(10.0)) + 1);
908 /* split fn and ext */
909 ext = std::strrchr(fn, '.');
912 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
915 /* insert e.g. '%03d' between fn and ext */
916 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
917 snew(fnout, std::strlen(buf)+1);
918 std::strcpy(fnout, buf);
923 static void ana_trans(t_clusters *clust, int nf,
924 const char *transfn, const char *ntransfn, FILE *log,
925 t_rgb rlo, t_rgb rhi, const gmx_output_env_t *oenv)
930 int i, ntranst, maxtrans;
933 snew(ntrans, clust->ncl);
934 snew(trans, clust->ncl);
935 snew(axis, clust->ncl);
936 for (i = 0; i < clust->ncl; i++)
939 snew(trans[i], clust->ncl);
943 for (i = 1; i < nf; i++)
945 if (clust->cl[i] != clust->cl[i-1])
948 ntrans[clust->cl[i-1]-1]++;
949 ntrans[clust->cl[i]-1]++;
950 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
951 maxtrans = static_cast<int>(std::max(static_cast<real>(maxtrans), trans[clust->cl[i]-1][clust->cl[i-1]-1]));
954 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
955 "max %d between two specific clusters\n", ntranst, maxtrans);
958 fp = gmx_ffopen(transfn, "w");
959 i = std::min(maxtrans+1, 80);
960 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
961 "from cluster", "to cluster",
962 clust->ncl, clust->ncl, axis, axis, trans,
963 0, maxtrans, rlo, rhi, &i);
968 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
970 for (i = 0; i < clust->ncl; i++)
972 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
977 for (i = 0; i < clust->ncl; i++)
985 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
986 int natom, t_atoms *atoms, rvec *xtps,
987 real *mass, rvec **xx, real *time,
988 int ifsize, int *fitidx,
989 int iosize, int *outidx,
990 const char *trxfn, const char *sizefn,
991 const char *transfn, const char *ntransfn,
992 const char *clustidfn, gmx_bool bAverage,
993 int write_ncl, int write_nst, real rmsmin,
994 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
995 const gmx_output_env_t *oenv)
997 FILE *size_fp = nullptr;
998 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
999 t_trxstatus *trxout = nullptr;
1000 t_trxstatus *trxsout = nullptr;
1001 int i, i1, cl, nstr, *structure, first = 0, midstr;
1002 gmx_bool *bWrite = nullptr;
1003 real r, clrmsd, midrmsd;
1004 rvec *xav = nullptr;
1009 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1013 /* do we write all structures? */
1016 trxsfn = parse_filename(trxfn, std::max(write_ncl, clust->ncl));
1019 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1020 bAverage ? "average" : "middle", trxfn);
1023 /* find out what we want to tell the user:
1024 Writing [all structures|structures with rmsd > %g] for
1025 {all|first %d} clusters {with more than %d structures} to %s */
1028 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1032 sprintf(buf1, "all structures");
1034 buf2[0] = buf3[0] = '\0';
1035 if (write_ncl >= clust->ncl)
1039 sprintf(buf2, "all ");
1044 sprintf(buf2, "the first %d ", write_ncl);
1048 sprintf(buf3, " with more than %d structures", write_nst);
1050 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1051 ffprintf(stderr, log, buf);
1054 /* Prepare a reference structure for the orientation of the clusters */
1057 reset_x(ifsize, fitidx, natom, nullptr, xtps, mass);
1059 trxout = open_trx(trxfn, "w");
1060 /* Calculate the average structure in each cluster, *
1061 * all structures are fitted to the first struture of the cluster */
1065 if (transfn || ntransfn)
1067 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1072 FILE *fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1073 if (output_env_get_print_xvgr_codes(oenv))
1075 fprintf(fp, "@ s0 symbol 2\n");
1076 fprintf(fp, "@ s0 symbol size 0.2\n");
1077 fprintf(fp, "@ s0 linestyle 0\n");
1079 for (i = 0; i < nf; i++)
1081 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1087 size_fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1088 if (output_env_get_print_xvgr_codes(oenv))
1090 fprintf(size_fp, "@g%d type %s\n", 0, "bar");
1093 snew(structure, nf);
1094 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1095 "cl.", "#st", "rmsd", "middle", "rmsd");
1096 for (cl = 1; cl <= clust->ncl; cl++)
1098 /* prepare structures (fit, middle, average) */
1101 for (i = 0; i < natom; i++)
1107 for (i1 = 0; i1 < nf; i1++)
1109 if (clust->cl[i1] == cl)
1111 structure[nstr] = i1;
1113 if (trxfn && (bAverage || write_ncl) )
1117 reset_x(ifsize, fitidx, natom, nullptr, xx[i1], mass);
1125 do_fit(natom, mass, xx[first], xx[i1]);
1129 for (i = 0; i < natom; i++)
1131 rvec_inc(xav[i], xx[i1][i]);
1139 fprintf(size_fp, "%8d %8d\n", cl, nstr);
1144 for (i1 = 0; i1 < nstr; i1++)
1149 for (i = 0; i < nstr; i++)
1153 r += rmsd[structure[i]][structure[i1]];
1157 r += rmsd[structure[i1]][structure[i]];
1164 midstr = structure[i1];
1171 /* dump cluster info to logfile */
1174 sprintf(buf1, "%6.3f", clrmsd);
1179 sprintf(buf2, "%5.3f", midrmsd);
1187 sprintf(buf1, "%5s", "");
1188 sprintf(buf2, "%5s", "");
1190 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1191 for (i = 0; i < nstr; i++)
1193 if ((i % 7 == 0) && i)
1195 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1202 fprintf(log, "%s %6g", buf, time[i1]);
1206 /* write structures to trajectory file(s) */
1211 for (i = 0; i < nstr; i++)
1216 if (cl < write_ncl+1 && nstr > write_nst)
1218 /* Dump all structures for this cluster */
1219 /* generate numbered filename (there is a %d in trxfn!) */
1220 sprintf(buf, trxsfn, cl);
1221 trxsout = open_trx(buf, "w");
1222 for (i = 0; i < nstr; i++)
1227 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1231 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1237 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1238 xx[structure[i]], nullptr, nullptr);
1243 /* Dump the average structure for this cluster */
1246 for (i = 0; i < natom; i++)
1248 svmul(1.0/nstr, xav[i], xav[i]);
1253 for (i = 0; i < natom; i++)
1255 copy_rvec(xx[midstr][i], xav[i]);
1259 reset_x(ifsize, fitidx, natom, nullptr, xav, mass);
1264 do_fit(natom, mass, xtps, xav);
1266 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, nullptr, nullptr);
1291 static void convert_mat(t_matrix *mat, t_mat *rms)
1296 matrix2real(mat, rms->mat);
1297 /* free input xpm matrix data */
1298 for (i = 0; i < mat->nx; i++)
1300 sfree(mat->matrix[i]);
1304 for (i = 0; i < mat->nx; i++)
1306 for (j = i; j < mat->nx; j++)
1308 rms->sumrms += rms->mat[i][j];
1309 rms->maxrms = std::max(rms->maxrms, rms->mat[i][j]);
1312 rms->minrms = std::min(rms->minrms, rms->mat[i][j]);
1319 int gmx_cluster(int argc, char *argv[])
1321 const char *desc[] = {
1322 "[THISMODULE] can cluster structures using several different methods.",
1323 "Distances between structures can be determined from a trajectory",
1324 "or read from an [REF].xpm[ref] matrix file with the [TT]-dm[tt] option.",
1325 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1326 "can be used to define the distance between structures.[PAR]",
1328 "single linkage: add a structure to a cluster when its distance to any",
1329 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1331 "Jarvis Patrick: add a structure to a cluster when this structure",
1332 "and a structure in the cluster have each other as neighbors and",
1333 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1334 "of a structure are the M closest structures or all structures within",
1335 "[TT]cutoff[tt].[PAR]",
1337 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1338 "the order of the frames is using the smallest possible increments.",
1339 "With this it is possible to make a smooth animation going from one",
1340 "structure to another with the largest possible (e.g.) RMSD between",
1341 "them, however the intermediate steps should be as small as possible.",
1342 "Applications could be to visualize a potential of mean force",
1343 "ensemble of simulations or a pulling simulation. Obviously the user",
1344 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1345 "The final result can be inspect visually by looking at the matrix",
1346 "[REF].xpm[ref] file, which should vary smoothly from bottom to top.[PAR]",
1348 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1350 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1351 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1352 "Count number of neighbors using cut-off, take structure with",
1353 "largest number of neighbors with all its neighbors as cluster",
1354 "and eliminate it from the pool of clusters. Repeat for remaining",
1355 "structures in pool.[PAR]",
1357 "When the clustering algorithm assigns each structure to exactly one",
1358 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1359 "file is supplied, the structure with",
1360 "the smallest average distance to the others or the average structure",
1361 "or all structures for each cluster will be written to a trajectory",
1362 "file. When writing all structures, separate numbered files are made",
1363 "for each cluster.[PAR]",
1365 "Two output files are always written:",
1367 " * [TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1368 " and a graphical depiction of the clusters in the lower right half",
1369 " When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1370 " when two structures are in the same cluster.",
1371 " When [TT]-minstruct[tt] > 1 different colors will be used for each",
1373 " * [TT]-g[tt] writes information on the options used and a detailed list",
1374 " of all clusters and their members.",
1377 "Additionally, a number of optional output files can be written:",
1379 " * [TT]-dist[tt] writes the RMSD distribution.",
1380 " * [TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1381 " diagonalization.",
1382 " * [TT]-sz[tt] writes the cluster sizes.",
1383 " * [TT]-tr[tt] writes a matrix of the number transitions between",
1385 " * [TT]-ntr[tt] writes the total number of transitions to or from",
1387 " * [TT]-clid[tt] writes the cluster number as a function of time.",
1388 " * [TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1389 " structure of each cluster or writes numbered files with cluster members",
1390 " for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1391 " [TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1392 " structure with the smallest average RMSD from all other structures",
1397 int nf, i, i1, i2, j;
1398 gmx_int64_t nrms = 0;
1401 rvec *xtps, *usextps, *x1, **xx = nullptr;
1402 const char *fn, *trx_out_fn;
1404 t_mat *rms, *orig = nullptr;
1409 t_matrix *readmat = nullptr;
1412 int isize = 0, ifsize = 0, iosize = 0;
1413 int *index = nullptr, *fitidx = nullptr, *outidx = nullptr;
1415 real rmsd, **d1, **d2, *time = nullptr, time_invfac, *mass = nullptr;
1416 char buf[STRLEN], buf1[80];
1417 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1419 int method, ncluster = 0;
1420 static const char *methodname[] = {
1421 nullptr, "linkage", "jarvis-patrick", "monte-carlo",
1422 "diagonalization", "gromos", nullptr
1425 m_null, m_linkage, m_jarvis_patrick,
1426 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1428 /* Set colors for plotting: white = zero RMS, black = maximum */
1429 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1430 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1431 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1432 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1433 static int nlevels = 40, skip = 1;
1434 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1435 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1436 static int niter = 10000, nrandom = 0, seed = 0, write_ncl = 0, write_nst = 1, minstruct = 1;
1437 static real kT = 1e-3;
1438 static int M = 10, P = 3;
1439 gmx_output_env_t *oenv;
1440 gmx_rmpbc_t gpbc = nullptr;
1443 { "-dista", FALSE, etBOOL, {&bRMSdist},
1444 "Use RMSD of distances instead of RMS deviation" },
1445 { "-nlevels", FALSE, etINT, {&nlevels},
1446 "Discretize RMSD matrix in this number of levels" },
1447 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1448 "RMSD cut-off (nm) for two structures to be neighbor" },
1449 { "-fit", FALSE, etBOOL, {&bFit},
1450 "Use least squares fitting before RMSD calculation" },
1451 { "-max", FALSE, etREAL, {&scalemax},
1452 "Maximum level in RMSD matrix" },
1453 { "-skip", FALSE, etINT, {&skip},
1454 "Only analyze every nr-th frame" },
1455 { "-av", FALSE, etBOOL, {&bAverage},
1456 "Write average instead of middle structure for each cluster" },
1457 { "-wcl", FALSE, etINT, {&write_ncl},
1458 "Write the structures for this number of clusters to numbered files" },
1459 { "-nst", FALSE, etINT, {&write_nst},
1460 "Only write all structures if more than this number of structures per cluster" },
1461 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1462 "minimum rms difference with rest of cluster for writing structures" },
1463 { "-method", FALSE, etENUM, {methodname},
1464 "Method for cluster determination" },
1465 { "-minstruct", FALSE, etINT, {&minstruct},
1466 "Minimum number of structures in cluster for coloring in the [REF].xpm[ref] file" },
1467 { "-binary", FALSE, etBOOL, {&bBinary},
1468 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1469 "is given by [TT]-cutoff[tt]" },
1470 { "-M", FALSE, etINT, {&M},
1471 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1472 "0 is use cutoff" },
1473 { "-P", FALSE, etINT, {&P},
1474 "Number of identical nearest neighbors required to form a cluster" },
1475 { "-seed", FALSE, etINT, {&seed},
1476 "Random number seed for Monte Carlo clustering algorithm (0 means generate)" },
1477 { "-niter", FALSE, etINT, {&niter},
1478 "Number of iterations for MC" },
1479 { "-nrandom", FALSE, etINT, {&nrandom},
1480 "The first iterations for MC may be done complete random, to shuffle the frames" },
1481 { "-kT", FALSE, etREAL, {&kT},
1482 "Boltzmann weighting factor for Monte Carlo optimization "
1483 "(zero turns off uphill steps)" },
1484 { "-pbc", FALSE, etBOOL,
1485 { &bPBC }, "PBC check" }
1488 { efTRX, "-f", nullptr, ffOPTRD },
1489 { efTPS, "-s", nullptr, ffOPTRD },
1490 { efNDX, nullptr, nullptr, ffOPTRD },
1491 { efXPM, "-dm", "rmsd", ffOPTRD },
1492 { efXPM, "-om", "rmsd-raw", ffWRITE },
1493 { efXPM, "-o", "rmsd-clust", ffWRITE },
1494 { efLOG, "-g", "cluster", ffWRITE },
1495 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1496 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1497 { efXVG, "-conv", "mc-conv", ffOPTWR },
1498 { efXVG, "-sz", "clust-size", ffOPTWR},
1499 { efXPM, "-tr", "clust-trans", ffOPTWR},
1500 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1501 { efXVG, "-clid", "clust-id", ffOPTWR},
1502 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1504 #define NFILE asize(fnm)
1506 if (!parse_common_args(&argc, argv,
1507 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
1508 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr,
1515 bReadMat = opt2bSet("-dm", NFILE, fnm);
1516 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1517 if (opt2parg_bSet("-av", asize(pa), pa) ||
1518 opt2parg_bSet("-wcl", asize(pa), pa) ||
1519 opt2parg_bSet("-nst", asize(pa), pa) ||
1520 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1521 opt2bSet("-cl", NFILE, fnm) )
1523 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1527 trx_out_fn = nullptr;
1529 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1532 "\nWarning: assuming the time unit in %s is %s\n",
1533 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv).c_str());
1535 if (trx_out_fn && !bReadTraj)
1537 fprintf(stderr, "\nWarning: "
1538 "cannot write cluster structures without reading trajectory\n"
1539 " ignoring option -cl %s\n", trx_out_fn);
1543 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1549 gmx_fatal(FARGS, "Invalid method");
1552 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1553 method == m_gromos );
1556 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1558 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1559 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1561 /* check input and write parameters to log file */
1562 bUseRmsdCut = FALSE;
1563 if (method == m_jarvis_patrick)
1565 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1566 if ((M < 0) || (M == 1))
1568 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1572 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1579 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1583 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1588 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1591 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1593 else /* method != m_jarvis */
1595 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1597 if (bUseRmsdCut && method != m_jarvis_patrick)
1599 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1601 if (method == m_monte_carlo)
1603 fprintf(log, "Using %d iterations\n", niter);
1608 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1614 /* don't read mass-database as masses (and top) are not used */
1615 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtps, nullptr, box,
1619 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1622 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1623 bReadMat ? "" : " and RMSD calculation");
1624 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1625 1, &ifsize, &fitidx, &grpname);
1628 fprintf(stderr, "\nSelect group for output:\n");
1629 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1630 1, &iosize, &outidx, &grpname);
1631 /* merge and convert both index groups: */
1632 /* first copy outidx to index. let outidx refer to elements in index */
1633 snew(index, iosize);
1635 for (i = 0; i < iosize; i++)
1637 index[i] = outidx[i];
1640 /* now lookup elements from fitidx in index, add them if necessary
1641 and also let fitidx refer to elements in index */
1642 for (i = 0; i < ifsize; i++)
1645 while (j < isize && index[j] != fitidx[i])
1651 /* slow this way, but doesn't matter much */
1653 srenew(index, isize);
1655 index[j] = fitidx[i];
1659 else /* !trx_out_fn */
1663 for (i = 0; i < ifsize; i++)
1665 index[i] = fitidx[i];
1673 /* Loop over first coordinate file */
1674 fn = opt2fn("-f", NFILE, fnm);
1676 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1677 output_env_conv_times(oenv, nf, time);
1678 if (!bRMSdist || bAnalyze)
1680 /* Center all frames on zero */
1682 for (i = 0; i < ifsize; i++)
1684 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1688 for (i = 0; i < nf; i++)
1690 reset_x(ifsize, fitidx, isize, nullptr, xx[i], mass);
1696 gmx_rmpbc_done(gpbc);
1702 fprintf(stderr, "Reading rms distance matrix ");
1703 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1704 fprintf(stderr, "\n");
1705 if (readmat[0].nx != readmat[0].ny)
1707 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1708 readmat[0].nx, readmat[0].ny);
1710 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1712 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1713 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1718 time = readmat[0].axis_x;
1719 time_invfac = output_env_get_time_invfactor(oenv);
1720 for (i = 0; i < nf; i++)
1722 time[i] *= time_invfac;
1725 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1726 convert_mat(&(readmat[0]), rms);
1728 nlevels = readmat[0].nmap;
1730 else /* !bReadMat */
1732 rms = init_mat(nf, method == m_diagonalize);
1733 nrms = (static_cast<gmx_int64_t>(nf)*static_cast<gmx_int64_t>(nf-1))/2;
1736 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1737 /* Initialize work array */
1739 for (i1 = 0; i1 < nf; i1++)
1741 for (i2 = i1+1; i2 < nf; i2++)
1743 for (i = 0; i < isize; i++)
1745 copy_rvec(xx[i1][i], x1[i]);
1749 do_fit(isize, mass, xx[i2], x1);
1751 rmsd = rmsdev(isize, mass, xx[i2], x1);
1752 set_mat_entry(rms, i1, i2, rmsd);
1755 fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
1762 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1764 /* Initiate work arrays */
1767 for (i = 0; (i < isize); i++)
1772 for (i1 = 0; i1 < nf; i1++)
1774 calc_dist(isize, xx[i1], d1);
1775 for (i2 = i1+1; (i2 < nf); i2++)
1777 calc_dist(isize, xx[i2], d2);
1778 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1781 fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
1784 /* Clean up work arrays */
1785 for (i = 0; (i < isize); i++)
1793 fprintf(stderr, "\n\n");
1795 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1796 rms->minrms, rms->maxrms);
1797 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1798 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1799 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1800 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1802 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1803 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1805 if (bAnalyze && (rmsmin < rms->minrms) )
1807 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1808 rmsmin, rms->minrms);
1810 if (bAnalyze && (rmsmin > rmsdcut) )
1812 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1816 /* Plot the rmsd distribution */
1817 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1821 for (i1 = 0; (i1 < nf); i1++)
1823 for (i2 = 0; (i2 < nf); i2++)
1825 if (rms->mat[i1][i2] < rmsdcut)
1827 rms->mat[i1][i2] = 0;
1831 rms->mat[i1][i2] = 1;
1841 /* Now sort the matrix and write it out again */
1842 gather(rms, rmsdcut, &clust);
1845 /* Do a diagonalization */
1846 snew(eigenvalues, nf);
1847 snew(eigenvectors, nf*nf);
1848 std::memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1849 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1850 sfree(eigenvectors);
1852 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1853 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1854 for (i = 0; (i < nf); i++)
1856 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1861 orig = init_mat(rms->nn, FALSE);
1863 copy_t_mat(orig, rms);
1864 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1865 opt2fn_null("-conv", NFILE, fnm), oenv);
1867 case m_jarvis_patrick:
1868 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1871 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1874 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1877 if (method == m_monte_carlo || method == m_diagonalize)
1879 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1887 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1891 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1893 init_t_atoms(&useatoms, isize, FALSE);
1894 snew(usextps, isize);
1895 useatoms.resinfo = top.atoms.resinfo;
1896 for (i = 0; i < isize; i++)
1898 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1899 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1900 useatoms.nres = std::max(useatoms.nres, useatoms.atom[i].resind+1);
1901 copy_rvec(xtps[index[i]], usextps[i]);
1903 useatoms.nr = isize;
1904 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1905 ifsize, fitidx, iosize, outidx,
1906 bReadTraj ? trx_out_fn : nullptr,
1907 opt2fn_null("-sz", NFILE, fnm),
1908 opt2fn_null("-tr", NFILE, fnm),
1909 opt2fn_null("-ntr", NFILE, fnm),
1910 opt2fn_null("-clid", NFILE, fnm),
1911 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1912 rlo_bot, rhi_bot, oenv);
1916 if (bBinary && !bAnalyze)
1918 /* Make the clustering visible */
1919 for (i2 = 0; (i2 < nf); i2++)
1921 for (i1 = i2+1; (i1 < nf); i1++)
1923 if (rms->mat[i1][i2])
1925 rms->mat[i1][i2] = rms->maxrms;
1931 fp = opt2FILE("-o", NFILE, fnm, "w");
1932 fprintf(stderr, "Writing rms distance/clustering matrix ");
1935 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1936 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1937 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1941 auto timeLabel = output_env_get_time_label(oenv);
1942 auto title = gmx::formatString("RMS%sDeviation / Cluster Index",
1943 bRMSdist ? " Distance " : " ");
1946 write_xpm_split(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1947 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1948 rlo_top, rhi_top, 0.0, ncluster,
1949 &ncluster, TRUE, rlo_bot, rhi_bot);
1953 write_xpm(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1954 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1955 rlo_top, rhi_top, &nlevels);
1958 fprintf(stderr, "\n");
1960 if (nullptr != orig)
1962 fp = opt2FILE("-om", NFILE, fnm, "w");
1963 auto timeLabel = output_env_get_time_label(oenv);
1964 auto title = gmx::formatString("RMS%sDeviation", bRMSdist ? " Distance " : " ");
1965 write_xpm(fp, 0, title, "RMSD (nm)", timeLabel, timeLabel,
1966 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1967 rlo_top, rhi_top, &nlevels);
1972 /* now show what we've done */
1973 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1974 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1975 if (method == m_diagonalize)
1977 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1979 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1982 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1983 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1984 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1986 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), nullptr);