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, 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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trnio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/cmat.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/legacyheaders/macros.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/legacyheaders/viewit.h"
54 #include "gromacs/linearalgebra/eigensolver.h"
55 #include "gromacs/math/do_fit.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/random/random.h"
60 #include "gromacs/topology/index.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/smalloc.h"
66 /* print to two file pointers at once (i.e. stderr and log) */
68 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
70 fprintf(fp1, "%s", buf);
71 fprintf(fp2, "%s", buf);
74 /* just print a prepared buffer to fp1 and fp2 */
76 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
78 lo_ffprintf(fp1, fp2, buf);
81 /* prepare buffer with one argument, then print to fp1 and fp2 */
83 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
85 sprintf(buf, fmt, arg);
86 lo_ffprintf(fp1, fp2, buf);
89 /* prepare buffer with one argument, then print to fp1 and fp2 */
91 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
93 sprintf(buf, fmt, arg);
94 lo_ffprintf(fp1, fp2, buf);
97 /* prepare buffer with one argument, then print to fp1 and fp2 */
99 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
101 sprintf(buf, fmt, arg);
102 lo_ffprintf(fp1, fp2, buf);
105 /* prepare buffer with two arguments, then print to fp1 and fp2 */
107 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
109 sprintf(buf, fmt, arg1, arg2);
110 lo_ffprintf(fp1, fp2, buf);
113 /* prepare buffer with two arguments, then print to fp1 and fp2 */
115 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
117 sprintf(buf, fmt, arg1, arg2);
118 lo_ffprintf(fp1, fp2, buf);
121 /* prepare buffer with two arguments, then print to fp1 and fp2 */
123 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
125 sprintf(buf, fmt, arg1, arg2);
126 lo_ffprintf(fp1, fp2, buf);
139 void pr_energy(FILE *fp, real e)
141 fprintf(fp, "Energy: %8.4f\n", e);
144 void cp_index(int nn, int from[], int to[])
148 for (i = 0; (i < nn); i++)
154 void mc_optimize(FILE *log, t_mat *m, real *time,
155 int maxiter, int nrandom,
157 const char *conv, output_env_t oenv)
160 real ecur, enext, emin, prob, enorm;
161 int i, j, iswap, jswap, nn, nuphill = 0;
167 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
170 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
171 printf("by reordering the frames to minimize the path between the two structures\n");
172 printf("that have the largest pairwise RMSD.\n");
175 enorm = m->mat[0][0];
176 for (i = 0; (i < m->n1); i++)
178 for (j = 0; (j < m->nn); j++)
180 if (m->mat[i][j] > enorm)
182 enorm = m->mat[i][j];
188 if ((iswap == -1) || (jswap == -1))
190 fprintf(stderr, "Matrix contains identical values in all fields\n");
193 swap_rows(m, 0, iswap);
194 swap_rows(m, m->n1-1, jswap);
195 emin = ecur = mat_energy(m);
196 printf("Largest distance %g between %d and %d. Energy: %g.\n",
197 enorm, iswap, jswap, emin);
199 rng = gmx_rng_init(seed);
202 /* Initiate and store global minimum */
203 minimum = init_mat(nn, m->b1D);
205 copy_t_mat(minimum, m);
209 fp = xvgropen(conv, "Convergence of the MC optimization",
210 "Energy", "Step", oenv);
212 for (i = 0; (i < maxiter); i++)
214 /* Generate new swapping candidates */
217 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
218 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
220 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
222 /* Apply swap and compute energy */
223 swap_rows(m, iswap, jswap);
224 enext = mat_energy(m);
226 /* Compute probability */
228 if ((enext < ecur) || (i < nrandom))
233 /* Store global minimum */
234 copy_t_mat(minimum, m);
240 /* Try Monte Carlo step */
241 prob = exp(-(enext-ecur)/(enorm*kT));
244 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
251 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
252 i, iswap, jswap, enext, prob);
255 fprintf(fp, "%6d %10g\n", i, enext);
261 swap_rows(m, jswap, iswap);
264 fprintf(log, "%d uphill steps were taken during optimization\n",
267 /* Now swap the matrix to get it into global minimum mode */
268 copy_t_mat(m, minimum);
270 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
271 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
272 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
273 for (i = 0; (i < m->nn); i++)
275 fprintf(log, "%10g %5d %10g\n",
278 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
287 static void calc_dist(int nind, rvec x[], real **d)
293 for (i = 0; (i < nind-1); i++)
296 for (j = i+1; (j < nind); j++)
298 /* Should use pbc_dx when analysing multiple molecueles,
299 * but the box is not stored for every frame.
301 rvec_sub(xi, x[j], dx);
307 static real rms_dist(int isize, real **d, real **d_r)
313 for (i = 0; (i < isize-1); i++)
315 for (j = i+1; (j < isize); j++)
317 r = d[i][j]-d_r[i][j];
321 r2 /= (isize*(isize-1))/2;
326 static int rms_dist_comp(const void *a, const void *b)
333 if (da->dist - db->dist < 0)
337 else if (da->dist - db->dist > 0)
344 static int clust_id_comp(const void *a, const void *b)
351 return da->clust - db->clust;
354 static int nrnb_comp(const void *a, const void *b)
361 /* return the b-a, we want highest first */
362 return db->nr - da->nr;
365 void gather(t_mat *m, real cutoff, t_clusters *clust)
369 int i, j, k, nn, cid, n1, diff;
372 /* First we sort the entries in the RMSD matrix */
376 for (i = k = 0; (i < n1); i++)
378 for (j = i+1; (j < n1); j++, k++)
382 d[k].dist = m->mat[i][j];
387 gmx_incons("gather algortihm");
389 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
391 /* Now we make a cluster index for all of the conformations */
394 /* Now we check the closest structures, and equalize their cluster numbers */
395 fprintf(stderr, "Linking structures ");
398 fprintf(stderr, "*");
400 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
402 diff = c[d[k].j].clust - c[d[k].i].clust;
408 c[d[k].j].clust = c[d[k].i].clust;
412 c[d[k].i].clust = c[d[k].j].clust;
418 fprintf(stderr, "\nSorting and renumbering clusters\n");
419 /* Sort on cluster number */
420 qsort(c, n1, sizeof(c[0]), clust_id_comp);
422 /* Renumber clusters */
424 for (k = 1; k < n1; k++)
426 if (c[k].clust != c[k-1].clust)
439 for (k = 0; (k < n1); k++)
441 fprintf(debug, "Cluster index for conformation %d: %d\n",
442 c[k].conf, c[k].clust);
446 for (k = 0; k < n1; k++)
448 clust->cl[c[k].conf] = c[k].clust;
455 gmx_bool jp_same(int **nnb, int i, int j, int P)
461 for (k = 0; nnb[i][k] >= 0; k++)
463 bIn = bIn || (nnb[i][k] == j);
471 for (k = 0; nnb[j][k] >= 0; k++)
473 bIn = bIn || (nnb[j][k] == i);
481 for (ii = 0; nnb[i][ii] >= 0; ii++)
483 for (jj = 0; nnb[j][jj] >= 0; jj++)
485 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
495 static void jarvis_patrick(int n1, real **mat, int M, int P,
496 real rmsdcut, t_clusters *clust)
501 int i, j, k, cid, diff, max;
510 /* First we sort the entries in the RMSD matrix row by row.
511 * This gives us the nearest neighbor list.
515 for (i = 0; (i < n1); i++)
517 for (j = 0; (j < n1); j++)
520 row[j].dist = mat[i][j];
522 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
525 /* Put the M nearest neighbors in the list */
527 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
531 nnb[i][k] = row[j].j;
539 /* Put all neighbors nearer than rmsdcut in the list */
542 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
551 nnb[i][k] = row[j].j;
557 srenew(nnb[i], max+1);
565 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
566 for (i = 0; (i < n1); i++)
568 fprintf(debug, "i:%5d nbs:", i);
569 for (j = 0; nnb[i][j] >= 0; j++)
571 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
573 fprintf(debug, "\n");
578 fprintf(stderr, "Linking structures ");
579 /* Use mcpy for temporary storage of booleans */
580 mcpy = mk_matrix(n1, n1, FALSE);
581 for (i = 0; i < n1; i++)
583 for (j = i+1; j < n1; j++)
585 mcpy[i][j] = jp_same(nnb, i, j, P);
590 fprintf(stderr, "*");
592 for (i = 0; i < n1; i++)
594 for (j = i+1; j < n1; j++)
598 diff = c[j].clust - c[i].clust;
604 c[j].clust = c[i].clust;
608 c[i].clust = c[j].clust;
617 fprintf(stderr, "\nSorting and renumbering clusters\n");
618 /* Sort on cluster number */
619 qsort(c, n1, sizeof(c[0]), clust_id_comp);
621 /* Renumber clusters */
623 for (k = 1; k < n1; k++)
625 if (c[k].clust != c[k-1].clust)
637 for (k = 0; k < n1; k++)
639 clust->cl[c[k].conf] = c[k].clust;
643 for (k = 0; (k < n1); k++)
645 fprintf(debug, "Cluster index for conformation %d: %d\n",
646 c[k].conf, c[k].clust);
650 /* Again, I don't see the point in this... (AF) */
651 /* for(i=0; (i<n1); i++) { */
652 /* for(j=0; (j<n1); j++) */
653 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
655 /* for(i=0; (i<n1); i++) { */
656 /* for(j=0; (j<n1); j++) */
657 /* mat[i][j] = mcpy[i][j]; */
659 done_matrix(n1, &mcpy);
662 for (i = 0; (i < n1); i++)
669 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
673 /* dump neighbor list */
674 fprintf(fp, "%s", title);
675 for (i = 0; (i < n1); i++)
677 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
678 for (j = 0; j < nnb[i].nr; j++)
680 fprintf(fp, "%5d", nnb[i].nb[j]);
686 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
690 int i, j, k, j1, max;
692 /* Put all neighbors nearer than rmsdcut in the list */
693 fprintf(stderr, "Making list of neighbors within cutoff ");
696 for (i = 0; (i < n1); i++)
700 /* put all neighbors within cut-off in list */
701 for (j = 0; j < n1; j++)
703 if (mat[i][j] < rmsdcut)
708 srenew(nnb[i].nb, max);
714 /* store nr of neighbors, we'll need that */
716 if (i%(1+n1/100) == 0)
718 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
721 fprintf(stderr, "%3d%%\n", 100);
724 /* sort neighbor list on number of neighbors, largest first */
725 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
729 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
732 /* turn first structure with all its neighbors (largest) into cluster
733 remove them from pool of structures and repeat for all remaining */
734 fprintf(stderr, "Finding clusters %4d", 0);
735 /* cluster id's start at 1: */
739 /* set cluster id (k) for first item in neighborlist */
740 for (j = 0; j < nnb[0].nr; j++)
742 clust->cl[nnb[0].nb[j]] = k;
748 /* adjust number of neighbors for others, taking removals into account: */
749 for (i = 1; i < n1 && nnb[i].nr; i++)
752 for (j = 0; j < nnb[i].nr; j++)
754 /* if this neighbor wasn't removed */
755 if (clust->cl[nnb[i].nb[j]] == 0)
757 /* shift the rest (j1<=j) */
758 nnb[i].nb[j1] = nnb[i].nb[j];
763 /* now j1 is the new number of neighbors */
766 /* sort again on nnb[].nr, because we have new # neighbors: */
767 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
768 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
770 fprintf(stderr, "\b\b\b\b%4d", k);
774 fprintf(stderr, "\n");
778 fprintf(debug, "Clusters (%d):\n", k);
779 for (i = 0; i < n1; i++)
781 fprintf(debug, " %3d", clust->cl[i]);
783 fprintf(debug, "\n");
789 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
790 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
795 int i, i0, j, max_nf;
803 natom = read_first_x(oenv, &status, fn, &t, &x, box);
810 gmx_rmpbc(gpbc, natom, box, x);
816 srenew(*time, max_nf);
821 /* Store only the interesting atoms */
822 for (j = 0; (j < isize); j++)
824 copy_rvec(x[index[j]], xx[i0][j]);
831 while (read_next_x(oenv, status, &t, x, box));
832 fprintf(stderr, "Allocated %lu bytes for frames\n",
833 (unsigned long) (max_nf*isize*sizeof(**xx)));
834 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
841 static int plot_clusters(int nf, real **mat, t_clusters *clust,
844 int i, j, ncluster, ci;
845 int *cl_id, *nstruct, *strind;
850 for (i = 0; i < nf; i++)
853 cl_id[i] = clust->cl[i];
857 for (i = 0; i < nf; i++)
859 if (nstruct[i] >= minstruct)
862 for (j = 0; (j < nf); j++)
866 strind[j] = ncluster;
872 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
873 ncluster, minstruct);
875 for (i = 0; (i < nf); i++)
878 for (j = 0; j < i; j++)
880 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
882 /* color different clusters with different colors, as long as
883 we don't run out of colors */
884 mat[i][j] = strind[i];
899 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
903 for (i = 0; i < nf; i++)
905 for (j = 0; j < i; j++)
907 if (clust->cl[i] == clust->cl[j])
919 static char *parse_filename(const char *fn, int maxnr)
928 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
930 /* number of digits needed in numbering */
931 i = (int)(log(maxnr)/log(10)) + 1;
932 /* split fn and ext */
933 ext = strrchr(fn, '.');
936 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
939 /* insert e.g. '%03d' between fn and ext */
940 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
941 snew(fnout, strlen(buf)+1);
947 static void ana_trans(t_clusters *clust, int nf,
948 const char *transfn, const char *ntransfn, FILE *log,
949 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
954 int i, ntranst, maxtrans;
957 snew(ntrans, clust->ncl);
958 snew(trans, clust->ncl);
959 snew(axis, clust->ncl);
960 for (i = 0; i < clust->ncl; i++)
963 snew(trans[i], clust->ncl);
967 for (i = 1; i < nf; i++)
969 if (clust->cl[i] != clust->cl[i-1])
972 ntrans[clust->cl[i-1]-1]++;
973 ntrans[clust->cl[i]-1]++;
974 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
975 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
978 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
979 "max %d between two specific clusters\n", ntranst, maxtrans);
982 fp = gmx_ffopen(transfn, "w");
983 i = min(maxtrans+1, 80);
984 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
985 "from cluster", "to cluster",
986 clust->ncl, clust->ncl, axis, axis, trans,
987 0, maxtrans, rlo, rhi, &i);
992 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
994 for (i = 0; i < clust->ncl; i++)
996 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1001 for (i = 0; i < clust->ncl; i++)
1009 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1010 int natom, t_atoms *atoms, rvec *xtps,
1011 real *mass, rvec **xx, real *time,
1012 int ifsize, atom_id *fitidx,
1013 int iosize, atom_id *outidx,
1014 const char *trxfn, const char *sizefn,
1015 const char *transfn, const char *ntransfn,
1016 const char *clustidfn, gmx_bool bAverage,
1017 int write_ncl, int write_nst, real rmsmin,
1018 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1019 const output_env_t oenv)
1022 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1023 t_trxstatus *trxout = NULL;
1024 t_trxstatus *trxsout = NULL;
1025 int i, i1, cl, nstr, *structure, first = 0, midstr;
1026 gmx_bool *bWrite = NULL;
1027 real r, clrmsd, midrmsd;
1033 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1037 /* do we write all structures? */
1040 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1043 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1044 bAverage ? "average" : "middle", trxfn);
1047 /* find out what we want to tell the user:
1048 Writing [all structures|structures with rmsd > %g] for
1049 {all|first %d} clusters {with more than %d structures} to %s */
1052 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1056 sprintf(buf1, "all structures");
1058 buf2[0] = buf3[0] = '\0';
1059 if (write_ncl >= clust->ncl)
1063 sprintf(buf2, "all ");
1068 sprintf(buf2, "the first %d ", write_ncl);
1072 sprintf(buf3, " with more than %d structures", write_nst);
1074 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1075 ffprintf(stderr, log, buf);
1078 /* Prepare a reference structure for the orientation of the clusters */
1081 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1083 trxout = open_trx(trxfn, "w");
1084 /* Calculate the average structure in each cluster, *
1085 * all structures are fitted to the first struture of the cluster */
1089 if (transfn || ntransfn)
1091 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1096 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1097 if (output_env_get_print_xvgr_codes(oenv))
1099 fprintf(fp, "@ s0 symbol 2\n");
1100 fprintf(fp, "@ s0 symbol size 0.2\n");
1101 fprintf(fp, "@ s0 linestyle 0\n");
1103 for (i = 0; i < nf; i++)
1105 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1111 /* FIXME: This file is never closed. */
1112 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1113 if (output_env_get_print_xvgr_codes(oenv))
1115 fprintf(fp, "@g%d type %s\n", 0, "bar");
1118 snew(structure, nf);
1119 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1120 "cl.", "#st", "rmsd", "middle", "rmsd");
1121 for (cl = 1; cl <= clust->ncl; cl++)
1123 /* prepare structures (fit, middle, average) */
1126 for (i = 0; i < natom; i++)
1132 for (i1 = 0; i1 < nf; i1++)
1134 if (clust->cl[i1] == cl)
1136 structure[nstr] = i1;
1138 if (trxfn && (bAverage || write_ncl) )
1142 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1150 do_fit(natom, mass, xx[first], xx[i1]);
1154 for (i = 0; i < natom; i++)
1156 rvec_inc(xav[i], xx[i1][i]);
1164 fprintf(fp, "%8d %8d\n", cl, nstr);
1169 for (i1 = 0; i1 < nstr; i1++)
1174 for (i = 0; i < nstr; i++)
1178 r += rmsd[structure[i]][structure[i1]];
1182 r += rmsd[structure[i1]][structure[i]];
1189 midstr = structure[i1];
1196 /* dump cluster info to logfile */
1199 sprintf(buf1, "%6.3f", clrmsd);
1204 sprintf(buf2, "%5.3f", midrmsd);
1212 sprintf(buf1, "%5s", "");
1213 sprintf(buf2, "%5s", "");
1215 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1216 for (i = 0; i < nstr; i++)
1218 if ((i % 7 == 0) && i)
1220 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1227 fprintf(log, "%s %6g", buf, time[i1]);
1231 /* write structures to trajectory file(s) */
1236 for (i = 0; i < nstr; i++)
1241 if (cl < write_ncl+1 && nstr > write_nst)
1243 /* Dump all structures for this cluster */
1244 /* generate numbered filename (there is a %d in trxfn!) */
1245 sprintf(buf, trxsfn, cl);
1246 trxsout = open_trx(buf, "w");
1247 for (i = 0; i < nstr; i++)
1252 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1256 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1262 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1263 xx[structure[i]], NULL, NULL);
1268 /* Dump the average structure for this cluster */
1271 for (i = 0; i < natom; i++)
1273 svmul(1.0/nstr, xav[i], xav[i]);
1278 for (i = 0; i < natom; i++)
1280 copy_rvec(xx[midstr][i], xav[i]);
1284 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1289 do_fit(natom, mass, xtps, xav);
1292 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1312 static void convert_mat(t_matrix *mat, t_mat *rms)
1317 matrix2real(mat, rms->mat);
1318 /* free input xpm matrix data */
1319 for (i = 0; i < mat->nx; i++)
1321 sfree(mat->matrix[i]);
1325 for (i = 0; i < mat->nx; i++)
1327 for (j = i; j < mat->nx; j++)
1329 rms->sumrms += rms->mat[i][j];
1330 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1333 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1340 int gmx_cluster(int argc, char *argv[])
1342 const char *desc[] = {
1343 "[THISMODULE] can cluster structures using several different methods.",
1344 "Distances between structures can be determined from a trajectory",
1345 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1346 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1347 "can be used to define the distance between structures.[PAR]",
1349 "single linkage: add a structure to a cluster when its distance to any",
1350 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1352 "Jarvis Patrick: add a structure to a cluster when this structure",
1353 "and a structure in the cluster have each other as neighbors and",
1354 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1355 "of a structure are the M closest structures or all structures within",
1356 "[TT]cutoff[tt].[PAR]",
1358 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1359 "the order of the frames is using the smallest possible increments.",
1360 "With this it is possible to make a smooth animation going from one",
1361 "structure to another with the largest possible (e.g.) RMSD between",
1362 "them, however the intermediate steps should be as small as possible.",
1363 "Applications could be to visualize a potential of mean force",
1364 "ensemble of simulations or a pulling simulation. Obviously the user",
1365 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1366 "The final result can be inspect visually by looking at the matrix",
1367 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1369 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1371 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1372 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1373 "Count number of neighbors using cut-off, take structure with",
1374 "largest number of neighbors with all its neighbors as cluster",
1375 "and eliminate it from the pool of clusters. Repeat for remaining",
1376 "structures in pool.[PAR]",
1378 "When the clustering algorithm assigns each structure to exactly one",
1379 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1380 "file is supplied, the structure with",
1381 "the smallest average distance to the others or the average structure",
1382 "or all structures for each cluster will be written to a trajectory",
1383 "file. When writing all structures, separate numbered files are made",
1384 "for each cluster.[PAR]",
1386 "Two output files are always written:[BR]",
1387 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1388 "and a graphical depiction of the clusters in the lower right half",
1389 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1390 "when two structures are in the same cluster.",
1391 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1393 "[TT]-g[tt] writes information on the options used and a detailed list",
1394 "of all clusters and their members.[PAR]",
1396 "Additionally, a number of optional output files can be written:[BR]",
1397 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1398 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1399 "diagonalization.[BR]",
1400 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1401 "[TT]-tr[tt] writes a matrix of the number transitions between",
1402 "cluster pairs.[BR]",
1403 "[TT]-ntr[tt] writes the total number of transitions to or from",
1404 "each cluster.[BR]",
1405 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1406 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1407 "structure of each cluster or writes numbered files with cluster members",
1408 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1409 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1410 "structure with the smallest average RMSD from all other structures",
1411 "of the cluster.[BR]",
1415 int nf, i, i1, i2, j;
1416 gmx_int64_t nrms = 0;
1419 rvec *xtps, *usextps, *x1, **xx = NULL;
1420 const char *fn, *trx_out_fn;
1422 t_mat *rms, *orig = NULL;
1427 t_matrix *readmat = NULL;
1430 int isize = 0, ifsize = 0, iosize = 0;
1431 atom_id *index = NULL, *fitidx, *outidx;
1433 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1434 char buf[STRLEN], buf1[80], title[STRLEN];
1435 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1437 int method, ncluster = 0;
1438 static const char *methodname[] = {
1439 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1440 "diagonalization", "gromos", NULL
1443 m_null, m_linkage, m_jarvis_patrick,
1444 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1446 /* Set colors for plotting: white = zero RMS, black = maximum */
1447 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1448 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1449 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1450 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1451 static int nlevels = 40, skip = 1;
1452 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1453 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1454 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1455 static real kT = 1e-3;
1456 static int M = 10, P = 3;
1458 gmx_rmpbc_t gpbc = NULL;
1461 { "-dista", FALSE, etBOOL, {&bRMSdist},
1462 "Use RMSD of distances instead of RMS deviation" },
1463 { "-nlevels", FALSE, etINT, {&nlevels},
1464 "Discretize RMSD matrix in this number of levels" },
1465 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1466 "RMSD cut-off (nm) for two structures to be neighbor" },
1467 { "-fit", FALSE, etBOOL, {&bFit},
1468 "Use least squares fitting before RMSD calculation" },
1469 { "-max", FALSE, etREAL, {&scalemax},
1470 "Maximum level in RMSD matrix" },
1471 { "-skip", FALSE, etINT, {&skip},
1472 "Only analyze every nr-th frame" },
1473 { "-av", FALSE, etBOOL, {&bAverage},
1474 "Write average iso middle structure for each cluster" },
1475 { "-wcl", FALSE, etINT, {&write_ncl},
1476 "Write the structures for this number of clusters to numbered files" },
1477 { "-nst", FALSE, etINT, {&write_nst},
1478 "Only write all structures if more than this number of structures per cluster" },
1479 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1480 "minimum rms difference with rest of cluster for writing structures" },
1481 { "-method", FALSE, etENUM, {methodname},
1482 "Method for cluster determination" },
1483 { "-minstruct", FALSE, etINT, {&minstruct},
1484 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1485 { "-binary", FALSE, etBOOL, {&bBinary},
1486 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1487 "is given by [TT]-cutoff[tt]" },
1488 { "-M", FALSE, etINT, {&M},
1489 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1490 "0 is use cutoff" },
1491 { "-P", FALSE, etINT, {&P},
1492 "Number of identical nearest neighbors required to form a cluster" },
1493 { "-seed", FALSE, etINT, {&seed},
1494 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1495 { "-niter", FALSE, etINT, {&niter},
1496 "Number of iterations for MC" },
1497 { "-nrandom", FALSE, etINT, {&nrandom},
1498 "The first iterations for MC may be done complete random, to shuffle the frames" },
1499 { "-kT", FALSE, etREAL, {&kT},
1500 "Boltzmann weighting factor for Monte Carlo optimization "
1501 "(zero turns off uphill steps)" },
1502 { "-pbc", FALSE, etBOOL,
1503 { &bPBC }, "PBC check" }
1506 { efTRX, "-f", NULL, ffOPTRD },
1507 { efTPS, "-s", NULL, ffOPTRD },
1508 { efNDX, NULL, NULL, ffOPTRD },
1509 { efXPM, "-dm", "rmsd", ffOPTRD },
1510 { efXPM, "-om", "rmsd-raw", ffWRITE },
1511 { efXPM, "-o", "rmsd-clust", ffWRITE },
1512 { efLOG, "-g", "cluster", ffWRITE },
1513 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1514 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1515 { efXVG, "-conv", "mc-conv", ffOPTWR },
1516 { efXVG, "-sz", "clust-size", ffOPTWR},
1517 { efXPM, "-tr", "clust-trans", ffOPTWR},
1518 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1519 { efXVG, "-clid", "clust-id", ffOPTWR},
1520 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1522 #define NFILE asize(fnm)
1524 if (!parse_common_args(&argc, argv,
1525 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
1526 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1533 bReadMat = opt2bSet("-dm", NFILE, fnm);
1534 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1535 if (opt2parg_bSet("-av", asize(pa), pa) ||
1536 opt2parg_bSet("-wcl", asize(pa), pa) ||
1537 opt2parg_bSet("-nst", asize(pa), pa) ||
1538 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1539 opt2bSet("-cl", NFILE, fnm) )
1541 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1547 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1550 "\nWarning: assuming the time unit in %s is %s\n",
1551 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1553 if (trx_out_fn && !bReadTraj)
1555 fprintf(stderr, "\nWarning: "
1556 "cannot write cluster structures without reading trajectory\n"
1557 " ignoring option -cl %s\n", trx_out_fn);
1561 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1567 gmx_fatal(FARGS, "Invalid method");
1570 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1571 method == m_gromos );
1574 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1576 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1577 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1579 /* check input and write parameters to log file */
1580 bUseRmsdCut = FALSE;
1581 if (method == m_jarvis_patrick)
1583 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1584 if ((M < 0) || (M == 1))
1586 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1590 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1597 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1601 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1606 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1609 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1611 else /* method != m_jarvis */
1613 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1615 if (bUseRmsdCut && method != m_jarvis_patrick)
1617 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1619 if (method == m_monte_carlo)
1621 fprintf(log, "Using %d iterations\n", niter);
1626 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1632 /* don't read mass-database as masses (and top) are not used */
1633 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1637 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1640 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1641 bReadMat ? "" : " and RMSD calculation");
1642 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1643 1, &ifsize, &fitidx, &grpname);
1646 fprintf(stderr, "\nSelect group for output:\n");
1647 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1648 1, &iosize, &outidx, &grpname);
1649 /* merge and convert both index groups: */
1650 /* first copy outidx to index. let outidx refer to elements in index */
1651 snew(index, iosize);
1653 for (i = 0; i < iosize; i++)
1655 index[i] = outidx[i];
1658 /* now lookup elements from fitidx in index, add them if necessary
1659 and also let fitidx refer to elements in index */
1660 for (i = 0; i < ifsize; i++)
1663 while (j < isize && index[j] != fitidx[i])
1669 /* slow this way, but doesn't matter much */
1671 srenew(index, isize);
1673 index[j] = fitidx[i];
1677 else /* !trx_out_fn */
1681 for (i = 0; i < ifsize; i++)
1683 index[i] = fitidx[i];
1691 /* Loop over first coordinate file */
1692 fn = opt2fn("-f", NFILE, fnm);
1694 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1695 output_env_conv_times(oenv, nf, time);
1696 if (!bRMSdist || bAnalyze)
1698 /* Center all frames on zero */
1700 for (i = 0; i < ifsize; i++)
1702 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1706 for (i = 0; i < nf; i++)
1708 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1714 gmx_rmpbc_done(gpbc);
1720 fprintf(stderr, "Reading rms distance matrix ");
1721 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1722 fprintf(stderr, "\n");
1723 if (readmat[0].nx != readmat[0].ny)
1725 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1726 readmat[0].nx, readmat[0].ny);
1728 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1730 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1731 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1736 time = readmat[0].axis_x;
1737 time_invfac = output_env_get_time_invfactor(oenv);
1738 for (i = 0; i < nf; i++)
1740 time[i] *= time_invfac;
1743 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1744 convert_mat(&(readmat[0]), rms);
1746 nlevels = readmat[0].nmap;
1748 else /* !bReadMat */
1750 rms = init_mat(nf, method == m_diagonalize);
1751 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1754 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1755 /* Initialize work array */
1757 for (i1 = 0; i1 < nf; i1++)
1759 for (i2 = i1+1; i2 < nf; i2++)
1761 for (i = 0; i < isize; i++)
1763 copy_rvec(xx[i1][i], x1[i]);
1767 do_fit(isize, mass, xx[i2], x1);
1769 rmsd = rmsdev(isize, mass, xx[i2], x1);
1770 set_mat_entry(rms, i1, i2, rmsd);
1772 nrms -= (gmx_int64_t) (nf-i1-1);
1773 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1779 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1781 /* Initiate work arrays */
1784 for (i = 0; (i < isize); i++)
1789 for (i1 = 0; i1 < nf; i1++)
1791 calc_dist(isize, xx[i1], d1);
1792 for (i2 = i1+1; (i2 < nf); i2++)
1794 calc_dist(isize, xx[i2], d2);
1795 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1798 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1800 /* Clean up work arrays */
1801 for (i = 0; (i < isize); i++)
1809 fprintf(stderr, "\n\n");
1811 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1812 rms->minrms, rms->maxrms);
1813 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1814 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1815 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1816 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1818 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1819 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1821 if (bAnalyze && (rmsmin < rms->minrms) )
1823 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1824 rmsmin, rms->minrms);
1826 if (bAnalyze && (rmsmin > rmsdcut) )
1828 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1832 /* Plot the rmsd distribution */
1833 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1837 for (i1 = 0; (i1 < nf); i1++)
1839 for (i2 = 0; (i2 < nf); i2++)
1841 if (rms->mat[i1][i2] < rmsdcut)
1843 rms->mat[i1][i2] = 0;
1847 rms->mat[i1][i2] = 1;
1857 /* Now sort the matrix and write it out again */
1858 gather(rms, rmsdcut, &clust);
1861 /* Do a diagonalization */
1862 snew(eigenvalues, nf);
1863 snew(eigenvectors, nf*nf);
1864 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1865 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1866 sfree(eigenvectors);
1868 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1869 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1870 for (i = 0; (i < nf); i++)
1872 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1877 orig = init_mat(rms->nn, FALSE);
1879 copy_t_mat(orig, rms);
1880 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1881 opt2fn_null("-conv", NFILE, fnm), oenv);
1883 case m_jarvis_patrick:
1884 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1887 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1890 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1893 if (method == m_monte_carlo || method == m_diagonalize)
1895 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1903 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1907 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1909 init_t_atoms(&useatoms, isize, FALSE);
1910 snew(usextps, isize);
1911 useatoms.resinfo = top.atoms.resinfo;
1912 for (i = 0; i < isize; i++)
1914 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1915 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1916 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1917 copy_rvec(xtps[index[i]], usextps[i]);
1919 useatoms.nr = isize;
1920 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1921 ifsize, fitidx, iosize, outidx,
1922 bReadTraj ? trx_out_fn : NULL,
1923 opt2fn_null("-sz", NFILE, fnm),
1924 opt2fn_null("-tr", NFILE, fnm),
1925 opt2fn_null("-ntr", NFILE, fnm),
1926 opt2fn_null("-clid", NFILE, fnm),
1927 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1928 rlo_bot, rhi_bot, oenv);
1932 if (bBinary && !bAnalyze)
1934 /* Make the clustering visible */
1935 for (i2 = 0; (i2 < nf); i2++)
1937 for (i1 = i2+1; (i1 < nf); i1++)
1939 if (rms->mat[i1][i2])
1941 rms->mat[i1][i2] = rms->maxrms;
1947 fp = opt2FILE("-o", NFILE, fnm, "w");
1948 fprintf(stderr, "Writing rms distance/clustering matrix ");
1951 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1952 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1953 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1957 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1958 sprintf(title, "RMS%sDeviation / Cluster Index",
1959 bRMSdist ? " Distance " : " ");
1962 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1963 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1964 rlo_top, rhi_top, 0.0, (real) ncluster,
1965 &ncluster, TRUE, rlo_bot, rhi_bot);
1969 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1970 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1971 rlo_top, rhi_top, &nlevels);
1974 fprintf(stderr, "\n");
1978 fp = opt2FILE("-om", NFILE, fnm, "w");
1979 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1980 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1981 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1982 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1983 rlo_top, rhi_top, &nlevels);
1988 /* now show what we've done */
1989 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1990 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1991 if (method == m_diagonalize)
1993 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1995 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1998 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1999 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
2000 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
2002 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);