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, 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/utility/smalloc.h"
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/tpxio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/random/random.h"
58 #include "gromacs/fileio/futil.h"
59 #include "gromacs/fileio/matio.h"
61 #include "gromacs/fileio/trnio.h"
65 #include "gromacs/linearalgebra/eigensolver.h"
66 #include "gromacs/math/do_fit.h"
67 #include "gromacs/legacyheaders/gmx_fatal.h"
69 /* print to two file pointers at once (i.e. stderr and log) */
71 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
73 fprintf(fp1, "%s", buf);
74 fprintf(fp2, "%s", buf);
77 /* just print a prepared buffer to fp1 and fp2 */
79 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
81 lo_ffprintf(fp1, fp2, buf);
84 /* prepare buffer with one argument, then print to fp1 and fp2 */
86 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
88 sprintf(buf, fmt, arg);
89 lo_ffprintf(fp1, fp2, buf);
92 /* prepare buffer with one argument, then print to fp1 and fp2 */
94 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
96 sprintf(buf, fmt, arg);
97 lo_ffprintf(fp1, fp2, buf);
100 /* prepare buffer with one argument, then print to fp1 and fp2 */
102 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
104 sprintf(buf, fmt, arg);
105 lo_ffprintf(fp1, fp2, buf);
108 /* prepare buffer with two arguments, then print to fp1 and fp2 */
110 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
112 sprintf(buf, fmt, arg1, arg2);
113 lo_ffprintf(fp1, fp2, buf);
116 /* prepare buffer with two arguments, then print to fp1 and fp2 */
118 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
120 sprintf(buf, fmt, arg1, arg2);
121 lo_ffprintf(fp1, fp2, buf);
124 /* prepare buffer with two arguments, then print to fp1 and fp2 */
126 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
128 sprintf(buf, fmt, arg1, arg2);
129 lo_ffprintf(fp1, fp2, buf);
142 void pr_energy(FILE *fp, real e)
144 fprintf(fp, "Energy: %8.4f\n", e);
147 void cp_index(int nn, int from[], int to[])
151 for (i = 0; (i < nn); i++)
157 void mc_optimize(FILE *log, t_mat *m, real *time,
158 int maxiter, int nrandom,
160 const char *conv, output_env_t oenv)
163 real ecur, enext, emin, prob, enorm;
164 int i, j, iswap, jswap, nn, nuphill = 0;
170 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
173 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
174 printf("by reordering the frames to minimize the path between the two structures\n");
175 printf("that have the largest pairwise RMSD.\n");
178 enorm = m->mat[0][0];
179 for (i = 0; (i < m->n1); i++)
181 for (j = 0; (j < m->nn); j++)
183 if (m->mat[i][j] > enorm)
185 enorm = m->mat[i][j];
191 if ((iswap == -1) || (jswap == -1))
193 fprintf(stderr, "Matrix contains identical values in all fields\n");
196 swap_rows(m, 0, iswap);
197 swap_rows(m, m->n1-1, jswap);
198 emin = ecur = mat_energy(m);
199 printf("Largest distance %g between %d and %d. Energy: %g.\n",
200 enorm, iswap, jswap, emin);
202 rng = gmx_rng_init(seed);
205 /* Initiate and store global minimum */
206 minimum = init_mat(nn, m->b1D);
208 copy_t_mat(minimum, m);
212 fp = xvgropen(conv, "Convergence of the MC optimization",
213 "Energy", "Step", oenv);
215 for (i = 0; (i < maxiter); i++)
217 /* Generate new swapping candidates */
220 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
221 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
223 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
225 /* Apply swap and compute energy */
226 swap_rows(m, iswap, jswap);
227 enext = mat_energy(m);
229 /* Compute probability */
231 if ((enext < ecur) || (i < nrandom))
236 /* Store global minimum */
237 copy_t_mat(minimum, m);
243 /* Try Monte Carlo step */
244 prob = exp(-(enext-ecur)/(enorm*kT));
247 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
254 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
255 i, iswap, jswap, enext, prob);
258 fprintf(fp, "%6d %10g\n", i, enext);
264 swap_rows(m, jswap, iswap);
267 fprintf(log, "%d uphill steps were taken during optimization\n",
270 /* Now swap the matrix to get it into global minimum mode */
271 copy_t_mat(m, minimum);
273 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
274 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
275 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
276 for (i = 0; (i < m->nn); i++)
278 fprintf(log, "%10g %5d %10g\n",
281 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
290 static void calc_dist(int nind, rvec x[], real **d)
296 for (i = 0; (i < nind-1); i++)
299 for (j = i+1; (j < nind); j++)
301 /* Should use pbc_dx when analysing multiple molecueles,
302 * but the box is not stored for every frame.
304 rvec_sub(xi, x[j], dx);
310 static real rms_dist(int isize, real **d, real **d_r)
316 for (i = 0; (i < isize-1); i++)
318 for (j = i+1; (j < isize); j++)
320 r = d[i][j]-d_r[i][j];
324 r2 /= (isize*(isize-1))/2;
329 static int rms_dist_comp(const void *a, const void *b)
336 if (da->dist - db->dist < 0)
340 else if (da->dist - db->dist > 0)
347 static int clust_id_comp(const void *a, const void *b)
354 return da->clust - db->clust;
357 static int nrnb_comp(const void *a, const void *b)
364 /* return the b-a, we want highest first */
365 return db->nr - da->nr;
368 void gather(t_mat *m, real cutoff, t_clusters *clust)
372 int i, j, k, nn, cid, n1, diff;
375 /* First we sort the entries in the RMSD matrix */
379 for (i = k = 0; (i < n1); i++)
381 for (j = i+1; (j < n1); j++, k++)
385 d[k].dist = m->mat[i][j];
390 gmx_incons("gather algortihm");
392 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
394 /* Now we make a cluster index for all of the conformations */
397 /* Now we check the closest structures, and equalize their cluster numbers */
398 fprintf(stderr, "Linking structures ");
401 fprintf(stderr, "*");
403 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
405 diff = c[d[k].j].clust - c[d[k].i].clust;
411 c[d[k].j].clust = c[d[k].i].clust;
415 c[d[k].i].clust = c[d[k].j].clust;
421 fprintf(stderr, "\nSorting and renumbering clusters\n");
422 /* Sort on cluster number */
423 qsort(c, n1, sizeof(c[0]), clust_id_comp);
425 /* Renumber clusters */
427 for (k = 1; k < n1; k++)
429 if (c[k].clust != c[k-1].clust)
442 for (k = 0; (k < n1); k++)
444 fprintf(debug, "Cluster index for conformation %d: %d\n",
445 c[k].conf, c[k].clust);
449 for (k = 0; k < n1; k++)
451 clust->cl[c[k].conf] = c[k].clust;
458 gmx_bool jp_same(int **nnb, int i, int j, int P)
464 for (k = 0; nnb[i][k] >= 0; k++)
466 bIn = bIn || (nnb[i][k] == j);
474 for (k = 0; nnb[j][k] >= 0; k++)
476 bIn = bIn || (nnb[j][k] == i);
484 for (ii = 0; nnb[i][ii] >= 0; ii++)
486 for (jj = 0; nnb[j][jj] >= 0; jj++)
488 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
498 static void jarvis_patrick(int n1, real **mat, int M, int P,
499 real rmsdcut, t_clusters *clust)
504 int i, j, k, cid, diff, max;
513 /* First we sort the entries in the RMSD matrix row by row.
514 * This gives us the nearest neighbor list.
518 for (i = 0; (i < n1); i++)
520 for (j = 0; (j < n1); j++)
523 row[j].dist = mat[i][j];
525 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
528 /* Put the M nearest neighbors in the list */
530 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
534 nnb[i][k] = row[j].j;
542 /* Put all neighbors nearer than rmsdcut in the list */
545 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
554 nnb[i][k] = row[j].j;
560 srenew(nnb[i], max+1);
568 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
569 for (i = 0; (i < n1); i++)
571 fprintf(debug, "i:%5d nbs:", i);
572 for (j = 0; nnb[i][j] >= 0; j++)
574 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
576 fprintf(debug, "\n");
581 fprintf(stderr, "Linking structures ");
582 /* Use mcpy for temporary storage of booleans */
583 mcpy = mk_matrix(n1, n1, FALSE);
584 for (i = 0; i < n1; i++)
586 for (j = i+1; j < n1; j++)
588 mcpy[i][j] = jp_same(nnb, i, j, P);
593 fprintf(stderr, "*");
595 for (i = 0; i < n1; i++)
597 for (j = i+1; j < n1; j++)
601 diff = c[j].clust - c[i].clust;
607 c[j].clust = c[i].clust;
611 c[i].clust = c[j].clust;
620 fprintf(stderr, "\nSorting and renumbering clusters\n");
621 /* Sort on cluster number */
622 qsort(c, n1, sizeof(c[0]), clust_id_comp);
624 /* Renumber clusters */
626 for (k = 1; k < n1; k++)
628 if (c[k].clust != c[k-1].clust)
640 for (k = 0; k < n1; k++)
642 clust->cl[c[k].conf] = c[k].clust;
646 for (k = 0; (k < n1); k++)
648 fprintf(debug, "Cluster index for conformation %d: %d\n",
649 c[k].conf, c[k].clust);
653 /* Again, I don't see the point in this... (AF) */
654 /* for(i=0; (i<n1); i++) { */
655 /* for(j=0; (j<n1); j++) */
656 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
658 /* for(i=0; (i<n1); i++) { */
659 /* for(j=0; (j<n1); j++) */
660 /* mat[i][j] = mcpy[i][j]; */
662 done_matrix(n1, &mcpy);
665 for (i = 0; (i < n1); i++)
672 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
676 /* dump neighbor list */
677 fprintf(fp, "%s", title);
678 for (i = 0; (i < n1); i++)
680 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
681 for (j = 0; j < nnb[i].nr; j++)
683 fprintf(fp, "%5d", nnb[i].nb[j]);
689 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
693 int i, j, k, j1, max;
695 /* Put all neighbors nearer than rmsdcut in the list */
696 fprintf(stderr, "Making list of neighbors within cutoff ");
699 for (i = 0; (i < n1); i++)
703 /* put all neighbors within cut-off in list */
704 for (j = 0; j < n1; j++)
706 if (mat[i][j] < rmsdcut)
711 srenew(nnb[i].nb, max);
717 /* store nr of neighbors, we'll need that */
719 if (i%(1+n1/100) == 0)
721 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
724 fprintf(stderr, "%3d%%\n", 100);
727 /* sort neighbor list on number of neighbors, largest first */
728 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
732 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
735 /* turn first structure with all its neighbors (largest) into cluster
736 remove them from pool of structures and repeat for all remaining */
737 fprintf(stderr, "Finding clusters %4d", 0);
738 /* cluster id's start at 1: */
742 /* set cluster id (k) for first item in neighborlist */
743 for (j = 0; j < nnb[0].nr; j++)
745 clust->cl[nnb[0].nb[j]] = k;
751 /* adjust number of neighbors for others, taking removals into account: */
752 for (i = 1; i < n1 && nnb[i].nr; i++)
755 for (j = 0; j < nnb[i].nr; j++)
757 /* if this neighbor wasn't removed */
758 if (clust->cl[nnb[i].nb[j]] == 0)
760 /* shift the rest (j1<=j) */
761 nnb[i].nb[j1] = nnb[i].nb[j];
766 /* now j1 is the new number of neighbors */
769 /* sort again on nnb[].nr, because we have new # neighbors: */
770 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
771 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
773 fprintf(stderr, "\b\b\b\b%4d", k);
777 fprintf(stderr, "\n");
781 fprintf(debug, "Clusters (%d):\n", k);
782 for (i = 0; i < n1; i++)
784 fprintf(debug, " %3d", clust->cl[i]);
786 fprintf(debug, "\n");
792 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
793 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
798 int i, i0, j, max_nf;
806 natom = read_first_x(oenv, &status, fn, &t, &x, box);
813 gmx_rmpbc(gpbc, natom, box, x);
819 srenew(*time, max_nf);
824 /* Store only the interesting atoms */
825 for (j = 0; (j < isize); j++)
827 copy_rvec(x[index[j]], xx[i0][j]);
834 while (read_next_x(oenv, status, &t, x, box));
835 fprintf(stderr, "Allocated %lu bytes for frames\n",
836 (unsigned long) (max_nf*isize*sizeof(**xx)));
837 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
844 static int plot_clusters(int nf, real **mat, t_clusters *clust,
847 int i, j, ncluster, ci;
848 int *cl_id, *nstruct, *strind;
853 for (i = 0; i < nf; i++)
856 cl_id[i] = clust->cl[i];
860 for (i = 0; i < nf; i++)
862 if (nstruct[i] >= minstruct)
865 for (j = 0; (j < nf); j++)
869 strind[j] = ncluster;
875 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
876 ncluster, minstruct);
878 for (i = 0; (i < nf); i++)
881 for (j = 0; j < i; j++)
883 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
885 /* color different clusters with different colors, as long as
886 we don't run out of colors */
887 mat[i][j] = strind[i];
902 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
906 for (i = 0; i < nf; i++)
908 for (j = 0; j < i; j++)
910 if (clust->cl[i] == clust->cl[j])
922 static char *parse_filename(const char *fn, int maxnr)
931 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
933 /* number of digits needed in numbering */
934 i = (int)(log(maxnr)/log(10)) + 1;
935 /* split fn and ext */
936 ext = strrchr(fn, '.');
939 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
942 /* insert e.g. '%03d' between fn and ext */
943 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
944 snew(fnout, strlen(buf)+1);
950 static void ana_trans(t_clusters *clust, int nf,
951 const char *transfn, const char *ntransfn, FILE *log,
952 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
957 int i, ntranst, maxtrans;
960 snew(ntrans, clust->ncl);
961 snew(trans, clust->ncl);
962 snew(axis, clust->ncl);
963 for (i = 0; i < clust->ncl; i++)
966 snew(trans[i], clust->ncl);
970 for (i = 1; i < nf; i++)
972 if (clust->cl[i] != clust->cl[i-1])
975 ntrans[clust->cl[i-1]-1]++;
976 ntrans[clust->cl[i]-1]++;
977 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
978 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
981 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
982 "max %d between two specific clusters\n", ntranst, maxtrans);
985 fp = gmx_ffopen(transfn, "w");
986 i = min(maxtrans+1, 80);
987 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
988 "from cluster", "to cluster",
989 clust->ncl, clust->ncl, axis, axis, trans,
990 0, maxtrans, rlo, rhi, &i);
995 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
997 for (i = 0; i < clust->ncl; i++)
999 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1004 for (i = 0; i < clust->ncl; i++)
1012 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1013 int natom, t_atoms *atoms, rvec *xtps,
1014 real *mass, rvec **xx, real *time,
1015 int ifsize, atom_id *fitidx,
1016 int iosize, atom_id *outidx,
1017 const char *trxfn, const char *sizefn,
1018 const char *transfn, const char *ntransfn,
1019 const char *clustidfn, gmx_bool bAverage,
1020 int write_ncl, int write_nst, real rmsmin,
1021 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1022 const output_env_t oenv)
1025 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1026 t_trxstatus *trxout = NULL;
1027 t_trxstatus *trxsout = NULL;
1028 int i, i1, cl, nstr, *structure, first = 0, midstr;
1029 gmx_bool *bWrite = NULL;
1030 real r, clrmsd, midrmsd;
1036 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1040 /* do we write all structures? */
1043 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1046 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1047 bAverage ? "average" : "middle", trxfn);
1050 /* find out what we want to tell the user:
1051 Writing [all structures|structures with rmsd > %g] for
1052 {all|first %d} clusters {with more than %d structures} to %s */
1055 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1059 sprintf(buf1, "all structures");
1061 buf2[0] = buf3[0] = '\0';
1062 if (write_ncl >= clust->ncl)
1066 sprintf(buf2, "all ");
1071 sprintf(buf2, "the first %d ", write_ncl);
1075 sprintf(buf3, " with more than %d structures", write_nst);
1077 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1078 ffprintf(stderr, log, buf);
1081 /* Prepare a reference structure for the orientation of the clusters */
1084 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1086 trxout = open_trx(trxfn, "w");
1087 /* Calculate the average structure in each cluster, *
1088 * all structures are fitted to the first struture of the cluster */
1092 if (transfn || ntransfn)
1094 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1099 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1100 fprintf(fp, "@ s0 symbol 2\n");
1101 fprintf(fp, "@ s0 symbol size 0.2\n");
1102 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 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1112 fprintf(fp, "@g%d type %s\n", 0, "bar");
1114 snew(structure, nf);
1115 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1116 "cl.", "#st", "rmsd", "middle", "rmsd");
1117 for (cl = 1; cl <= clust->ncl; cl++)
1119 /* prepare structures (fit, middle, average) */
1122 for (i = 0; i < natom; i++)
1128 for (i1 = 0; i1 < nf; i1++)
1130 if (clust->cl[i1] == cl)
1132 structure[nstr] = i1;
1134 if (trxfn && (bAverage || write_ncl) )
1138 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1146 do_fit(natom, mass, xx[first], xx[i1]);
1150 for (i = 0; i < natom; i++)
1152 rvec_inc(xav[i], xx[i1][i]);
1160 fprintf(fp, "%8d %8d\n", cl, nstr);
1165 for (i1 = 0; i1 < nstr; i1++)
1170 for (i = 0; i < nstr; i++)
1174 r += rmsd[structure[i]][structure[i1]];
1178 r += rmsd[structure[i1]][structure[i]];
1185 midstr = structure[i1];
1192 /* dump cluster info to logfile */
1195 sprintf(buf1, "%6.3f", clrmsd);
1200 sprintf(buf2, "%5.3f", midrmsd);
1208 sprintf(buf1, "%5s", "");
1209 sprintf(buf2, "%5s", "");
1211 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1212 for (i = 0; i < nstr; i++)
1214 if ((i % 7 == 0) && i)
1216 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1223 fprintf(log, "%s %6g", buf, time[i1]);
1227 /* write structures to trajectory file(s) */
1232 for (i = 0; i < nstr; i++)
1237 if (cl < write_ncl+1 && nstr > write_nst)
1239 /* Dump all structures for this cluster */
1240 /* generate numbered filename (there is a %d in trxfn!) */
1241 sprintf(buf, trxsfn, cl);
1242 trxsout = open_trx(buf, "w");
1243 for (i = 0; i < nstr; i++)
1248 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1252 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1258 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1259 xx[structure[i]], NULL, NULL);
1264 /* Dump the average structure for this cluster */
1267 for (i = 0; i < natom; i++)
1269 svmul(1.0/nstr, xav[i], xav[i]);
1274 for (i = 0; i < natom; i++)
1276 copy_rvec(xx[midstr][i], xav[i]);
1280 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1285 do_fit(natom, mass, xtps, xav);
1288 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1308 static void convert_mat(t_matrix *mat, t_mat *rms)
1313 matrix2real(mat, rms->mat);
1314 /* free input xpm matrix data */
1315 for (i = 0; i < mat->nx; i++)
1317 sfree(mat->matrix[i]);
1321 for (i = 0; i < mat->nx; i++)
1323 for (j = i; j < mat->nx; j++)
1325 rms->sumrms += rms->mat[i][j];
1326 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1329 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1336 int gmx_cluster(int argc, char *argv[])
1338 const char *desc[] = {
1339 "[THISMODULE] can cluster structures using several different methods.",
1340 "Distances between structures can be determined from a trajectory",
1341 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1342 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1343 "can be used to define the distance between structures.[PAR]",
1345 "single linkage: add a structure to a cluster when its distance to any",
1346 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1348 "Jarvis Patrick: add a structure to a cluster when this structure",
1349 "and a structure in the cluster have each other as neighbors and",
1350 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1351 "of a structure are the M closest structures or all structures within",
1352 "[TT]cutoff[tt].[PAR]",
1354 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1355 "the order of the frames is using the smallest possible increments.",
1356 "With this it is possible to make a smooth animation going from one",
1357 "structure to another with the largest possible (e.g.) RMSD between",
1358 "them, however the intermediate steps should be as small as possible.",
1359 "Applications could be to visualize a potential of mean force",
1360 "ensemble of simulations or a pulling simulation. Obviously the user",
1361 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1362 "The final result can be inspect visually by looking at the matrix",
1363 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1365 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1367 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1368 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1369 "Count number of neighbors using cut-off, take structure with",
1370 "largest number of neighbors with all its neighbors as cluster",
1371 "and eliminate it from the pool of clusters. Repeat for remaining",
1372 "structures in pool.[PAR]",
1374 "When the clustering algorithm assigns each structure to exactly one",
1375 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1376 "file is supplied, the structure with",
1377 "the smallest average distance to the others or the average structure",
1378 "or all structures for each cluster will be written to a trajectory",
1379 "file. When writing all structures, separate numbered files are made",
1380 "for each cluster.[PAR]",
1382 "Two output files are always written:[BR]",
1383 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1384 "and a graphical depiction of the clusters in the lower right half",
1385 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1386 "when two structures are in the same cluster.",
1387 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1389 "[TT]-g[tt] writes information on the options used and a detailed list",
1390 "of all clusters and their members.[PAR]",
1392 "Additionally, a number of optional output files can be written:[BR]",
1393 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1394 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1395 "diagonalization.[BR]",
1396 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1397 "[TT]-tr[tt] writes a matrix of the number transitions between",
1398 "cluster pairs.[BR]",
1399 "[TT]-ntr[tt] writes the total number of transitions to or from",
1400 "each cluster.[BR]",
1401 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1402 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1403 "structure of each cluster or writes numbered files with cluster members",
1404 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1405 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1406 "structure with the smallest average RMSD from all other structures",
1407 "of the cluster.[BR]",
1411 int nf, i, i1, i2, j;
1412 gmx_int64_t nrms = 0;
1415 rvec *xtps, *usextps, *x1, **xx = NULL;
1416 const char *fn, *trx_out_fn;
1418 t_mat *rms, *orig = NULL;
1423 t_matrix *readmat = NULL;
1426 int isize = 0, ifsize = 0, iosize = 0;
1427 atom_id *index = NULL, *fitidx, *outidx;
1429 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1430 char buf[STRLEN], buf1[80], title[STRLEN];
1431 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1433 int method, ncluster = 0;
1434 static const char *methodname[] = {
1435 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1436 "diagonalization", "gromos", NULL
1439 m_null, m_linkage, m_jarvis_patrick,
1440 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1442 /* Set colors for plotting: white = zero RMS, black = maximum */
1443 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1444 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1445 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1446 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1447 static int nlevels = 40, skip = 1;
1448 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1449 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1450 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1451 static real kT = 1e-3;
1452 static int M = 10, P = 3;
1454 gmx_rmpbc_t gpbc = NULL;
1457 { "-dista", FALSE, etBOOL, {&bRMSdist},
1458 "Use RMSD of distances instead of RMS deviation" },
1459 { "-nlevels", FALSE, etINT, {&nlevels},
1460 "Discretize RMSD matrix in this number of levels" },
1461 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1462 "RMSD cut-off (nm) for two structures to be neighbor" },
1463 { "-fit", FALSE, etBOOL, {&bFit},
1464 "Use least squares fitting before RMSD calculation" },
1465 { "-max", FALSE, etREAL, {&scalemax},
1466 "Maximum level in RMSD matrix" },
1467 { "-skip", FALSE, etINT, {&skip},
1468 "Only analyze every nr-th frame" },
1469 { "-av", FALSE, etBOOL, {&bAverage},
1470 "Write average iso middle structure for each cluster" },
1471 { "-wcl", FALSE, etINT, {&write_ncl},
1472 "Write the structures for this number of clusters to numbered files" },
1473 { "-nst", FALSE, etINT, {&write_nst},
1474 "Only write all structures if more than this number of structures per cluster" },
1475 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1476 "minimum rms difference with rest of cluster for writing structures" },
1477 { "-method", FALSE, etENUM, {methodname},
1478 "Method for cluster determination" },
1479 { "-minstruct", FALSE, etINT, {&minstruct},
1480 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1481 { "-binary", FALSE, etBOOL, {&bBinary},
1482 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1483 "is given by [TT]-cutoff[tt]" },
1484 { "-M", FALSE, etINT, {&M},
1485 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1486 "0 is use cutoff" },
1487 { "-P", FALSE, etINT, {&P},
1488 "Number of identical nearest neighbors required to form a cluster" },
1489 { "-seed", FALSE, etINT, {&seed},
1490 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1491 { "-niter", FALSE, etINT, {&niter},
1492 "Number of iterations for MC" },
1493 { "-nrandom", FALSE, etINT, {&nrandom},
1494 "The first iterations for MC may be done complete random, to shuffle the frames" },
1495 { "-kT", FALSE, etREAL, {&kT},
1496 "Boltzmann weighting factor for Monte Carlo optimization "
1497 "(zero turns off uphill steps)" },
1498 { "-pbc", FALSE, etBOOL,
1499 { &bPBC }, "PBC check" }
1502 { efTRX, "-f", NULL, ffOPTRD },
1503 { efTPS, "-s", NULL, ffOPTRD },
1504 { efNDX, NULL, NULL, ffOPTRD },
1505 { efXPM, "-dm", "rmsd", ffOPTRD },
1506 { efXPM, "-om", "rmsd-raw", ffWRITE },
1507 { efXPM, "-o", "rmsd-clust", ffWRITE },
1508 { efLOG, "-g", "cluster", ffWRITE },
1509 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1510 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1511 { efXVG, "-conv", "mc-conv", ffOPTWR },
1512 { efXVG, "-sz", "clust-size", ffOPTWR},
1513 { efXPM, "-tr", "clust-trans", ffOPTWR},
1514 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1515 { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1516 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1518 #define NFILE asize(fnm)
1520 if (!parse_common_args(&argc, argv,
1521 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1522 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1529 bReadMat = opt2bSet("-dm", NFILE, fnm);
1530 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1531 if (opt2parg_bSet("-av", asize(pa), pa) ||
1532 opt2parg_bSet("-wcl", asize(pa), pa) ||
1533 opt2parg_bSet("-nst", asize(pa), pa) ||
1534 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1535 opt2bSet("-cl", NFILE, fnm) )
1537 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1543 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1546 "\nWarning: assuming the time unit in %s is %s\n",
1547 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1549 if (trx_out_fn && !bReadTraj)
1551 fprintf(stderr, "\nWarning: "
1552 "cannot write cluster structures without reading trajectory\n"
1553 " ignoring option -cl %s\n", trx_out_fn);
1557 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1563 gmx_fatal(FARGS, "Invalid method");
1566 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1567 method == m_gromos );
1570 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1572 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1573 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1575 /* check input and write parameters to log file */
1576 bUseRmsdCut = FALSE;
1577 if (method == m_jarvis_patrick)
1579 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1580 if ((M < 0) || (M == 1))
1582 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1586 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1593 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1597 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1602 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1605 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1607 else /* method != m_jarvis */
1609 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1611 if (bUseRmsdCut && method != m_jarvis_patrick)
1613 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1615 if (method == m_monte_carlo)
1617 fprintf(log, "Using %d iterations\n", niter);
1622 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1628 /* don't read mass-database as masses (and top) are not used */
1629 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1633 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1636 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1637 bReadMat ? "" : " and RMSD calculation");
1638 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1639 1, &ifsize, &fitidx, &grpname);
1642 fprintf(stderr, "\nSelect group for output:\n");
1643 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1644 1, &iosize, &outidx, &grpname);
1645 /* merge and convert both index groups: */
1646 /* first copy outidx to index. let outidx refer to elements in index */
1647 snew(index, iosize);
1649 for (i = 0; i < iosize; i++)
1651 index[i] = outidx[i];
1654 /* now lookup elements from fitidx in index, add them if necessary
1655 and also let fitidx refer to elements in index */
1656 for (i = 0; i < ifsize; i++)
1659 while (j < isize && index[j] != fitidx[i])
1665 /* slow this way, but doesn't matter much */
1667 srenew(index, isize);
1669 index[j] = fitidx[i];
1673 else /* !trx_out_fn */
1677 for (i = 0; i < ifsize; i++)
1679 index[i] = fitidx[i];
1687 /* Loop over first coordinate file */
1688 fn = opt2fn("-f", NFILE, fnm);
1690 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1691 output_env_conv_times(oenv, nf, time);
1692 if (!bRMSdist || bAnalyze)
1694 /* Center all frames on zero */
1696 for (i = 0; i < ifsize; i++)
1698 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1702 for (i = 0; i < nf; i++)
1704 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1710 gmx_rmpbc_done(gpbc);
1716 fprintf(stderr, "Reading rms distance matrix ");
1717 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1718 fprintf(stderr, "\n");
1719 if (readmat[0].nx != readmat[0].ny)
1721 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1722 readmat[0].nx, readmat[0].ny);
1724 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1726 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1727 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1732 time = readmat[0].axis_x;
1733 time_invfac = output_env_get_time_invfactor(oenv);
1734 for (i = 0; i < nf; i++)
1736 time[i] *= time_invfac;
1739 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1740 convert_mat(&(readmat[0]), rms);
1742 nlevels = readmat[0].nmap;
1744 else /* !bReadMat */
1746 rms = init_mat(nf, method == m_diagonalize);
1747 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1750 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1751 /* Initialize work array */
1753 for (i1 = 0; i1 < nf; i1++)
1755 for (i2 = i1+1; i2 < nf; i2++)
1757 for (i = 0; i < isize; i++)
1759 copy_rvec(xx[i1][i], x1[i]);
1763 do_fit(isize, mass, xx[i2], x1);
1765 rmsd = rmsdev(isize, mass, xx[i2], x1);
1766 set_mat_entry(rms, i1, i2, rmsd);
1768 nrms -= (gmx_int64_t) (nf-i1-1);
1769 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1775 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1777 /* Initiate work arrays */
1780 for (i = 0; (i < isize); i++)
1785 for (i1 = 0; i1 < nf; i1++)
1787 calc_dist(isize, xx[i1], d1);
1788 for (i2 = i1+1; (i2 < nf); i2++)
1790 calc_dist(isize, xx[i2], d2);
1791 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1794 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1796 /* Clean up work arrays */
1797 for (i = 0; (i < isize); i++)
1805 fprintf(stderr, "\n\n");
1807 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1808 rms->minrms, rms->maxrms);
1809 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1810 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1811 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1812 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1814 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1815 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1817 if (bAnalyze && (rmsmin < rms->minrms) )
1819 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1820 rmsmin, rms->minrms);
1822 if (bAnalyze && (rmsmin > rmsdcut) )
1824 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1828 /* Plot the rmsd distribution */
1829 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1833 for (i1 = 0; (i1 < nf); i1++)
1835 for (i2 = 0; (i2 < nf); i2++)
1837 if (rms->mat[i1][i2] < rmsdcut)
1839 rms->mat[i1][i2] = 0;
1843 rms->mat[i1][i2] = 1;
1853 /* Now sort the matrix and write it out again */
1854 gather(rms, rmsdcut, &clust);
1857 /* Do a diagonalization */
1858 snew(eigenvalues, nf);
1859 snew(eigenvectors, nf*nf);
1860 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1861 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1862 sfree(eigenvectors);
1864 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1865 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1866 for (i = 0; (i < nf); i++)
1868 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1873 orig = init_mat(rms->nn, FALSE);
1875 copy_t_mat(orig, rms);
1876 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1877 opt2fn_null("-conv", NFILE, fnm), oenv);
1879 case m_jarvis_patrick:
1880 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1883 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1886 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1889 if (method == m_monte_carlo || method == m_diagonalize)
1891 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1899 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1903 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1905 init_t_atoms(&useatoms, isize, FALSE);
1906 snew(usextps, isize);
1907 useatoms.resinfo = top.atoms.resinfo;
1908 for (i = 0; i < isize; i++)
1910 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1911 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1912 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1913 copy_rvec(xtps[index[i]], usextps[i]);
1915 useatoms.nr = isize;
1916 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1917 ifsize, fitidx, iosize, outidx,
1918 bReadTraj ? trx_out_fn : NULL,
1919 opt2fn_null("-sz", NFILE, fnm),
1920 opt2fn_null("-tr", NFILE, fnm),
1921 opt2fn_null("-ntr", NFILE, fnm),
1922 opt2fn_null("-clid", NFILE, fnm),
1923 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1924 rlo_bot, rhi_bot, oenv);
1928 if (bBinary && !bAnalyze)
1930 /* Make the clustering visible */
1931 for (i2 = 0; (i2 < nf); i2++)
1933 for (i1 = i2+1; (i1 < nf); i1++)
1935 if (rms->mat[i1][i2])
1937 rms->mat[i1][i2] = rms->maxrms;
1943 fp = opt2FILE("-o", NFILE, fnm, "w");
1944 fprintf(stderr, "Writing rms distance/clustering matrix ");
1947 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1948 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1949 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1953 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1954 sprintf(title, "RMS%sDeviation / Cluster Index",
1955 bRMSdist ? " Distance " : " ");
1958 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1959 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1960 rlo_top, rhi_top, 0.0, (real) ncluster,
1961 &ncluster, TRUE, rlo_bot, rhi_bot);
1965 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1966 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1967 rlo_top, rhi_top, &nlevels);
1970 fprintf(stderr, "\n");
1974 fp = opt2FILE("-om", NFILE, fnm, "w");
1975 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1976 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1977 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1978 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1979 rlo_top, rhi_top, &nlevels);
1984 /* now show what we've done */
1985 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1986 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1987 if (method == m_diagonalize)
1989 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1991 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1994 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1995 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1996 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1998 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);