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.
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/fileio/tpxio.h"
48 #include "gromacs/fileio/trxio.h"
53 #include "gromacs/random/random.h"
57 #include "gromacs/fileio/futil.h"
58 #include "gromacs/fileio/matio.h"
61 #include "gromacs/fileio/trnio.h"
65 #include "gromacs/linearalgebra/eigensolver.h"
67 /* print to two file pointers at once (i.e. stderr and log) */
69 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
71 fprintf(fp1, "%s", buf);
72 fprintf(fp2, "%s", buf);
75 /* just print a prepared buffer to fp1 and fp2 */
77 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
79 lo_ffprintf(fp1, fp2, buf);
82 /* prepare buffer with one argument, then print to fp1 and fp2 */
84 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
86 sprintf(buf, fmt, arg);
87 lo_ffprintf(fp1, fp2, buf);
90 /* prepare buffer with one argument, then print to fp1 and fp2 */
92 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
94 sprintf(buf, fmt, arg);
95 lo_ffprintf(fp1, fp2, buf);
98 /* prepare buffer with one argument, then print to fp1 and fp2 */
100 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
102 sprintf(buf, fmt, arg);
103 lo_ffprintf(fp1, fp2, buf);
106 /* prepare buffer with two arguments, then print to fp1 and fp2 */
108 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
110 sprintf(buf, fmt, arg1, arg2);
111 lo_ffprintf(fp1, fp2, buf);
114 /* prepare buffer with two arguments, then print to fp1 and fp2 */
116 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
118 sprintf(buf, fmt, arg1, arg2);
119 lo_ffprintf(fp1, fp2, buf);
122 /* prepare buffer with two arguments, then print to fp1 and fp2 */
124 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
126 sprintf(buf, fmt, arg1, arg2);
127 lo_ffprintf(fp1, fp2, buf);
140 void pr_energy(FILE *fp, real e)
142 fprintf(fp, "Energy: %8.4f\n", e);
145 void cp_index(int nn, int from[], int to[])
149 for (i = 0; (i < nn); i++)
155 void mc_optimize(FILE *log, t_mat *m, real *time,
156 int maxiter, int nrandom,
158 const char *conv, output_env_t oenv)
161 real ecur, enext, emin, prob, enorm;
162 int i, j, iswap, jswap, nn, nuphill = 0;
168 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
171 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
172 printf("by reordering the frames to minimize the path between the two structures\n");
173 printf("that have the largest pairwise RMSD.\n");
176 enorm = m->mat[0][0];
177 for (i = 0; (i < m->n1); i++)
179 for (j = 0; (j < m->nn); j++)
181 if (m->mat[i][j] > enorm)
183 enorm = m->mat[i][j];
189 if ((iswap == -1) || (jswap == -1))
191 fprintf(stderr, "Matrix contains identical values in all fields\n");
194 swap_rows(m, 0, iswap);
195 swap_rows(m, m->n1-1, jswap);
196 emin = ecur = mat_energy(m);
197 printf("Largest distance %g between %d and %d. Energy: %g.\n",
198 enorm, iswap, jswap, emin);
200 rng = gmx_rng_init(seed);
203 /* Initiate and store global minimum */
204 minimum = init_mat(nn, m->b1D);
206 copy_t_mat(minimum, m);
210 fp = xvgropen(conv, "Convergence of the MC optimization",
211 "Energy", "Step", oenv);
213 for (i = 0; (i < maxiter); i++)
215 /* Generate new swapping candidates */
218 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
219 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
221 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
223 /* Apply swap and compute energy */
224 swap_rows(m, iswap, jswap);
225 enext = mat_energy(m);
227 /* Compute probability */
229 if ((enext < ecur) || (i < nrandom))
234 /* Store global minimum */
235 copy_t_mat(minimum, m);
241 /* Try Monte Carlo step */
242 prob = exp(-(enext-ecur)/(enorm*kT));
245 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
252 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
253 i, iswap, jswap, enext, prob);
256 fprintf(fp, "%6d %10g\n", i, enext);
262 swap_rows(m, jswap, iswap);
265 fprintf(log, "%d uphill steps were taken during optimization\n",
268 /* Now swap the matrix to get it into global minimum mode */
269 copy_t_mat(m, minimum);
271 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
272 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
273 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
274 for (i = 0; (i < m->nn); i++)
276 fprintf(log, "%10g %5d %10g\n",
279 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
288 static void calc_dist(int nind, rvec x[], real **d)
294 for (i = 0; (i < nind-1); i++)
297 for (j = i+1; (j < nind); j++)
299 /* Should use pbc_dx when analysing multiple molecueles,
300 * but the box is not stored for every frame.
302 rvec_sub(xi, x[j], dx);
308 static real rms_dist(int isize, real **d, real **d_r)
314 for (i = 0; (i < isize-1); i++)
316 for (j = i+1; (j < isize); j++)
318 r = d[i][j]-d_r[i][j];
322 r2 /= (isize*(isize-1))/2;
327 static int rms_dist_comp(const void *a, const void *b)
334 if (da->dist - db->dist < 0)
338 else if (da->dist - db->dist > 0)
345 static int clust_id_comp(const void *a, const void *b)
352 return da->clust - db->clust;
355 static int nrnb_comp(const void *a, const void *b)
362 /* return the b-a, we want highest first */
363 return db->nr - da->nr;
366 void gather(t_mat *m, real cutoff, t_clusters *clust)
370 int i, j, k, nn, cid, n1, diff;
373 /* First we sort the entries in the RMSD matrix */
377 for (i = k = 0; (i < n1); i++)
379 for (j = i+1; (j < n1); j++, k++)
383 d[k].dist = m->mat[i][j];
388 gmx_incons("gather algortihm");
390 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
392 /* Now we make a cluster index for all of the conformations */
395 /* Now we check the closest structures, and equalize their cluster numbers */
396 fprintf(stderr, "Linking structures ");
399 fprintf(stderr, "*");
401 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
403 diff = c[d[k].j].clust - c[d[k].i].clust;
409 c[d[k].j].clust = c[d[k].i].clust;
413 c[d[k].i].clust = c[d[k].j].clust;
419 fprintf(stderr, "\nSorting and renumbering clusters\n");
420 /* Sort on cluster number */
421 qsort(c, n1, sizeof(c[0]), clust_id_comp);
423 /* Renumber clusters */
425 for (k = 1; k < n1; k++)
427 if (c[k].clust != c[k-1].clust)
440 for (k = 0; (k < n1); k++)
442 fprintf(debug, "Cluster index for conformation %d: %d\n",
443 c[k].conf, c[k].clust);
447 for (k = 0; k < n1; k++)
449 clust->cl[c[k].conf] = c[k].clust;
456 gmx_bool jp_same(int **nnb, int i, int j, int P)
462 for (k = 0; nnb[i][k] >= 0; k++)
464 bIn = bIn || (nnb[i][k] == j);
472 for (k = 0; nnb[j][k] >= 0; k++)
474 bIn = bIn || (nnb[j][k] == i);
482 for (ii = 0; nnb[i][ii] >= 0; ii++)
484 for (jj = 0; nnb[j][jj] >= 0; jj++)
486 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
496 static void jarvis_patrick(int n1, real **mat, int M, int P,
497 real rmsdcut, t_clusters *clust)
502 int i, j, k, cid, diff, max;
511 /* First we sort the entries in the RMSD matrix row by row.
512 * This gives us the nearest neighbor list.
516 for (i = 0; (i < n1); i++)
518 for (j = 0; (j < n1); j++)
521 row[j].dist = mat[i][j];
523 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
526 /* Put the M nearest neighbors in the list */
528 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
532 nnb[i][k] = row[j].j;
540 /* Put all neighbors nearer than rmsdcut in the list */
543 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
552 nnb[i][k] = row[j].j;
558 srenew(nnb[i], max+1);
566 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
567 for (i = 0; (i < n1); i++)
569 fprintf(debug, "i:%5d nbs:", i);
570 for (j = 0; nnb[i][j] >= 0; j++)
572 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
574 fprintf(debug, "\n");
579 fprintf(stderr, "Linking structures ");
580 /* Use mcpy for temporary storage of booleans */
581 mcpy = mk_matrix(n1, n1, FALSE);
582 for (i = 0; i < n1; i++)
584 for (j = i+1; j < n1; j++)
586 mcpy[i][j] = jp_same(nnb, i, j, P);
591 fprintf(stderr, "*");
593 for (i = 0; i < n1; i++)
595 for (j = i+1; j < n1; j++)
599 diff = c[j].clust - c[i].clust;
605 c[j].clust = c[i].clust;
609 c[i].clust = c[j].clust;
618 fprintf(stderr, "\nSorting and renumbering clusters\n");
619 /* Sort on cluster number */
620 qsort(c, n1, sizeof(c[0]), clust_id_comp);
622 /* Renumber clusters */
624 for (k = 1; k < n1; k++)
626 if (c[k].clust != c[k-1].clust)
638 for (k = 0; k < n1; k++)
640 clust->cl[c[k].conf] = c[k].clust;
644 for (k = 0; (k < n1); k++)
646 fprintf(debug, "Cluster index for conformation %d: %d\n",
647 c[k].conf, c[k].clust);
651 /* Again, I don't see the point in this... (AF) */
652 /* for(i=0; (i<n1); i++) { */
653 /* for(j=0; (j<n1); j++) */
654 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
656 /* for(i=0; (i<n1); i++) { */
657 /* for(j=0; (j<n1); j++) */
658 /* mat[i][j] = mcpy[i][j]; */
660 done_matrix(n1, &mcpy);
663 for (i = 0; (i < n1); i++)
670 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
674 /* dump neighbor list */
675 fprintf(fp, "%s", title);
676 for (i = 0; (i < n1); i++)
678 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
679 for (j = 0; j < nnb[i].nr; j++)
681 fprintf(fp, "%5d", nnb[i].nb[j]);
687 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
691 int i, j, k, j1, max;
693 /* Put all neighbors nearer than rmsdcut in the list */
694 fprintf(stderr, "Making list of neighbors within cutoff ");
697 for (i = 0; (i < n1); i++)
701 /* put all neighbors within cut-off in list */
702 for (j = 0; j < n1; j++)
704 if (mat[i][j] < rmsdcut)
709 srenew(nnb[i].nb, max);
715 /* store nr of neighbors, we'll need that */
717 if (i%(1+n1/100) == 0)
719 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
722 fprintf(stderr, "%3d%%\n", 100);
725 /* sort neighbor list on number of neighbors, largest first */
726 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
730 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
733 /* turn first structure with all its neighbors (largest) into cluster
734 remove them from pool of structures and repeat for all remaining */
735 fprintf(stderr, "Finding clusters %4d", 0);
736 /* cluster id's start at 1: */
740 /* set cluster id (k) for first item in neighborlist */
741 for (j = 0; j < nnb[0].nr; j++)
743 clust->cl[nnb[0].nb[j]] = k;
749 /* adjust number of neighbors for others, taking removals into account: */
750 for (i = 1; i < n1 && nnb[i].nr; i++)
753 for (j = 0; j < nnb[i].nr; j++)
755 /* if this neighbor wasn't removed */
756 if (clust->cl[nnb[i].nb[j]] == 0)
758 /* shift the rest (j1<=j) */
759 nnb[i].nb[j1] = nnb[i].nb[j];
764 /* now j1 is the new number of neighbors */
767 /* sort again on nnb[].nr, because we have new # neighbors: */
768 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
769 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
771 fprintf(stderr, "\b\b\b\b%4d", k);
775 fprintf(stderr, "\n");
779 fprintf(debug, "Clusters (%d):\n", k);
780 for (i = 0; i < n1; i++)
782 fprintf(debug, " %3d", clust->cl[i]);
784 fprintf(debug, "\n");
790 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
791 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
796 int i, i0, j, max_nf;
804 natom = read_first_x(oenv, &status, fn, &t, &x, box);
811 gmx_rmpbc(gpbc, natom, box, x);
817 srenew(*time, max_nf);
822 /* Store only the interesting atoms */
823 for (j = 0; (j < isize); j++)
825 copy_rvec(x[index[j]], xx[i0][j]);
832 while (read_next_x(oenv, status, &t, x, box));
833 fprintf(stderr, "Allocated %lu bytes for frames\n",
834 (unsigned long) (max_nf*isize*sizeof(**xx)));
835 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
842 static int plot_clusters(int nf, real **mat, t_clusters *clust,
845 int i, j, ncluster, ci;
846 int *cl_id, *nstruct, *strind;
851 for (i = 0; i < nf; i++)
854 cl_id[i] = clust->cl[i];
858 for (i = 0; i < nf; i++)
860 if (nstruct[i] >= minstruct)
863 for (j = 0; (j < nf); j++)
867 strind[j] = ncluster;
873 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
874 ncluster, minstruct);
876 for (i = 0; (i < nf); i++)
879 for (j = 0; j < i; j++)
881 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
883 /* color different clusters with different colors, as long as
884 we don't run out of colors */
885 mat[i][j] = strind[i];
900 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
904 for (i = 0; i < nf; i++)
906 for (j = 0; j < i; j++)
908 if (clust->cl[i] == clust->cl[j])
920 static char *parse_filename(const char *fn, int maxnr)
929 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
931 /* number of digits needed in numbering */
932 i = (int)(log(maxnr)/log(10)) + 1;
933 /* split fn and ext */
934 ext = strrchr(fn, '.');
937 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
940 /* insert e.g. '%03d' between fn and ext */
941 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
942 snew(fnout, strlen(buf)+1);
948 static void ana_trans(t_clusters *clust, int nf,
949 const char *transfn, const char *ntransfn, FILE *log,
950 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
955 int i, ntranst, maxtrans;
958 snew(ntrans, clust->ncl);
959 snew(trans, clust->ncl);
960 snew(axis, clust->ncl);
961 for (i = 0; i < clust->ncl; i++)
964 snew(trans[i], clust->ncl);
968 for (i = 1; i < nf; i++)
970 if (clust->cl[i] != clust->cl[i-1])
973 ntrans[clust->cl[i-1]-1]++;
974 ntrans[clust->cl[i]-1]++;
975 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
976 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
979 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
980 "max %d between two specific clusters\n", ntranst, maxtrans);
983 fp = gmx_ffopen(transfn, "w");
984 i = min(maxtrans+1, 80);
985 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
986 "from cluster", "to cluster",
987 clust->ncl, clust->ncl, axis, axis, trans,
988 0, maxtrans, rlo, rhi, &i);
993 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
995 for (i = 0; i < clust->ncl; i++)
997 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1002 for (i = 0; i < clust->ncl; i++)
1010 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1011 int natom, t_atoms *atoms, rvec *xtps,
1012 real *mass, rvec **xx, real *time,
1013 int ifsize, atom_id *fitidx,
1014 int iosize, atom_id *outidx,
1015 const char *trxfn, const char *sizefn,
1016 const char *transfn, const char *ntransfn,
1017 const char *clustidfn, gmx_bool bAverage,
1018 int write_ncl, int write_nst, real rmsmin,
1019 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1020 const output_env_t oenv)
1023 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1024 t_trxstatus *trxout = NULL;
1025 t_trxstatus *trxsout = NULL;
1026 int i, i1, cl, nstr, *structure, first = 0, midstr;
1027 gmx_bool *bWrite = NULL;
1028 real r, clrmsd, midrmsd;
1034 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1038 /* do we write all structures? */
1041 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1044 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1045 bAverage ? "average" : "middle", trxfn);
1048 /* find out what we want to tell the user:
1049 Writing [all structures|structures with rmsd > %g] for
1050 {all|first %d} clusters {with more than %d structures} to %s */
1053 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1057 sprintf(buf1, "all structures");
1059 buf2[0] = buf3[0] = '\0';
1060 if (write_ncl >= clust->ncl)
1064 sprintf(buf2, "all ");
1069 sprintf(buf2, "the first %d ", write_ncl);
1073 sprintf(buf3, " with more than %d structures", write_nst);
1075 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1076 ffprintf(stderr, log, buf);
1079 /* Prepare a reference structure for the orientation of the clusters */
1082 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1084 trxout = open_trx(trxfn, "w");
1085 /* Calculate the average structure in each cluster, *
1086 * all structures are fitted to the first struture of the cluster */
1090 if (transfn || ntransfn)
1092 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1097 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1098 fprintf(fp, "@ s0 symbol 2\n");
1099 fprintf(fp, "@ s0 symbol size 0.2\n");
1100 fprintf(fp, "@ s0 linestyle 0\n");
1101 for (i = 0; i < nf; i++)
1103 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1109 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1110 fprintf(fp, "@g%d type %s\n", 0, "bar");
1112 snew(structure, nf);
1113 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1114 "cl.", "#st", "rmsd", "middle", "rmsd");
1115 for (cl = 1; cl <= clust->ncl; cl++)
1117 /* prepare structures (fit, middle, average) */
1120 for (i = 0; i < natom; i++)
1126 for (i1 = 0; i1 < nf; i1++)
1128 if (clust->cl[i1] == cl)
1130 structure[nstr] = i1;
1132 if (trxfn && (bAverage || write_ncl) )
1136 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1144 do_fit(natom, mass, xx[first], xx[i1]);
1148 for (i = 0; i < natom; i++)
1150 rvec_inc(xav[i], xx[i1][i]);
1158 fprintf(fp, "%8d %8d\n", cl, nstr);
1163 for (i1 = 0; i1 < nstr; i1++)
1168 for (i = 0; i < nstr; i++)
1172 r += rmsd[structure[i]][structure[i1]];
1176 r += rmsd[structure[i1]][structure[i]];
1183 midstr = structure[i1];
1190 /* dump cluster info to logfile */
1193 sprintf(buf1, "%6.3f", clrmsd);
1198 sprintf(buf2, "%5.3f", midrmsd);
1206 sprintf(buf1, "%5s", "");
1207 sprintf(buf2, "%5s", "");
1209 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1210 for (i = 0; i < nstr; i++)
1212 if ((i % 7 == 0) && i)
1214 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1221 fprintf(log, "%s %6g", buf, time[i1]);
1225 /* write structures to trajectory file(s) */
1230 for (i = 0; i < nstr; i++)
1235 if (cl < write_ncl+1 && nstr > write_nst)
1237 /* Dump all structures for this cluster */
1238 /* generate numbered filename (there is a %d in trxfn!) */
1239 sprintf(buf, trxsfn, cl);
1240 trxsout = open_trx(buf, "w");
1241 for (i = 0; i < nstr; i++)
1246 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1250 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1256 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1257 xx[structure[i]], NULL, NULL);
1262 /* Dump the average structure for this cluster */
1265 for (i = 0; i < natom; i++)
1267 svmul(1.0/nstr, xav[i], xav[i]);
1272 for (i = 0; i < natom; i++)
1274 copy_rvec(xx[midstr][i], xav[i]);
1278 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1283 do_fit(natom, mass, xtps, xav);
1286 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1306 static void convert_mat(t_matrix *mat, t_mat *rms)
1311 matrix2real(mat, rms->mat);
1312 /* free input xpm matrix data */
1313 for (i = 0; i < mat->nx; i++)
1315 sfree(mat->matrix[i]);
1319 for (i = 0; i < mat->nx; i++)
1321 for (j = i; j < mat->nx; j++)
1323 rms->sumrms += rms->mat[i][j];
1324 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1327 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1334 int gmx_cluster(int argc, char *argv[])
1336 const char *desc[] = {
1337 "[THISMODULE] can cluster structures using several different methods.",
1338 "Distances between structures can be determined from a trajectory",
1339 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1340 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1341 "can be used to define the distance between structures.[PAR]",
1343 "single linkage: add a structure to a cluster when its distance to any",
1344 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1346 "Jarvis Patrick: add a structure to a cluster when this structure",
1347 "and a structure in the cluster have each other as neighbors and",
1348 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1349 "of a structure are the M closest structures or all structures within",
1350 "[TT]cutoff[tt].[PAR]",
1352 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1353 "the order of the frames is using the smallest possible increments.",
1354 "With this it is possible to make a smooth animation going from one",
1355 "structure to another with the largest possible (e.g.) RMSD between",
1356 "them, however the intermediate steps should be as small as possible.",
1357 "Applications could be to visualize a potential of mean force",
1358 "ensemble of simulations or a pulling simulation. Obviously the user",
1359 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1360 "The final result can be inspect visually by looking at the matrix",
1361 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1363 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1365 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1366 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1367 "Count number of neighbors using cut-off, take structure with",
1368 "largest number of neighbors with all its neighbors as cluster",
1369 "and eliminate it from the pool of clusters. Repeat for remaining",
1370 "structures in pool.[PAR]",
1372 "When the clustering algorithm assigns each structure to exactly one",
1373 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1374 "file is supplied, the structure with",
1375 "the smallest average distance to the others or the average structure",
1376 "or all structures for each cluster will be written to a trajectory",
1377 "file. When writing all structures, separate numbered files are made",
1378 "for each cluster.[PAR]",
1380 "Two output files are always written:[BR]",
1381 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1382 "and a graphical depiction of the clusters in the lower right half",
1383 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1384 "when two structures are in the same cluster.",
1385 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1387 "[TT]-g[tt] writes information on the options used and a detailed list",
1388 "of all clusters and their members.[PAR]",
1390 "Additionally, a number of optional output files can be written:[BR]",
1391 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1392 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1393 "diagonalization.[BR]",
1394 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1395 "[TT]-tr[tt] writes a matrix of the number transitions between",
1396 "cluster pairs.[BR]",
1397 "[TT]-ntr[tt] writes the total number of transitions to or from",
1398 "each cluster.[BR]",
1399 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1400 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1401 "structure of each cluster or writes numbered files with cluster members",
1402 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1403 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1404 "structure with the smallest average RMSD from all other structures",
1405 "of the cluster.[BR]",
1409 int nf, i, i1, i2, j;
1410 gmx_int64_t nrms = 0;
1413 rvec *xtps, *usextps, *x1, **xx = NULL;
1414 const char *fn, *trx_out_fn;
1416 t_mat *rms, *orig = NULL;
1421 t_matrix *readmat = NULL;
1424 int isize = 0, ifsize = 0, iosize = 0;
1425 atom_id *index = NULL, *fitidx, *outidx;
1427 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1428 char buf[STRLEN], buf1[80], title[STRLEN];
1429 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1431 int method, ncluster = 0;
1432 static const char *methodname[] = {
1433 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1434 "diagonalization", "gromos", NULL
1437 m_null, m_linkage, m_jarvis_patrick,
1438 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1440 /* Set colors for plotting: white = zero RMS, black = maximum */
1441 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1442 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1443 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1444 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1445 static int nlevels = 40, skip = 1;
1446 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1447 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1448 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1449 static real kT = 1e-3;
1450 static int M = 10, P = 3;
1452 gmx_rmpbc_t gpbc = NULL;
1455 { "-dista", FALSE, etBOOL, {&bRMSdist},
1456 "Use RMSD of distances instead of RMS deviation" },
1457 { "-nlevels", FALSE, etINT, {&nlevels},
1458 "Discretize RMSD matrix in this number of levels" },
1459 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1460 "RMSD cut-off (nm) for two structures to be neighbor" },
1461 { "-fit", FALSE, etBOOL, {&bFit},
1462 "Use least squares fitting before RMSD calculation" },
1463 { "-max", FALSE, etREAL, {&scalemax},
1464 "Maximum level in RMSD matrix" },
1465 { "-skip", FALSE, etINT, {&skip},
1466 "Only analyze every nr-th frame" },
1467 { "-av", FALSE, etBOOL, {&bAverage},
1468 "Write average iso middle structure for each cluster" },
1469 { "-wcl", FALSE, etINT, {&write_ncl},
1470 "Write the structures for this number of clusters to numbered files" },
1471 { "-nst", FALSE, etINT, {&write_nst},
1472 "Only write all structures if more than this number of structures per cluster" },
1473 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1474 "minimum rms difference with rest of cluster for writing structures" },
1475 { "-method", FALSE, etENUM, {methodname},
1476 "Method for cluster determination" },
1477 { "-minstruct", FALSE, etINT, {&minstruct},
1478 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1479 { "-binary", FALSE, etBOOL, {&bBinary},
1480 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1481 "is given by [TT]-cutoff[tt]" },
1482 { "-M", FALSE, etINT, {&M},
1483 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1484 "0 is use cutoff" },
1485 { "-P", FALSE, etINT, {&P},
1486 "Number of identical nearest neighbors required to form a cluster" },
1487 { "-seed", FALSE, etINT, {&seed},
1488 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1489 { "-niter", FALSE, etINT, {&niter},
1490 "Number of iterations for MC" },
1491 { "-nrandom", FALSE, etINT, {&nrandom},
1492 "The first iterations for MC may be done complete random, to shuffle the frames" },
1493 { "-kT", FALSE, etREAL, {&kT},
1494 "Boltzmann weighting factor for Monte Carlo optimization "
1495 "(zero turns off uphill steps)" },
1496 { "-pbc", FALSE, etBOOL,
1497 { &bPBC }, "PBC check" }
1500 { efTRX, "-f", NULL, ffOPTRD },
1501 { efTPS, "-s", NULL, ffOPTRD },
1502 { efNDX, NULL, NULL, ffOPTRD },
1503 { efXPM, "-dm", "rmsd", ffOPTRD },
1504 { efXPM, "-om", "rmsd-raw", ffWRITE },
1505 { efXPM, "-o", "rmsd-clust", ffWRITE },
1506 { efLOG, "-g", "cluster", ffWRITE },
1507 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1508 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1509 { efXVG, "-conv", "mc-conv", ffOPTWR },
1510 { efXVG, "-sz", "clust-size", ffOPTWR},
1511 { efXPM, "-tr", "clust-trans", ffOPTWR},
1512 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1513 { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1514 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1516 #define NFILE asize(fnm)
1518 if (!parse_common_args(&argc, argv,
1519 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1520 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1527 bReadMat = opt2bSet("-dm", NFILE, fnm);
1528 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1529 if (opt2parg_bSet("-av", asize(pa), pa) ||
1530 opt2parg_bSet("-wcl", asize(pa), pa) ||
1531 opt2parg_bSet("-nst", asize(pa), pa) ||
1532 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1533 opt2bSet("-cl", NFILE, fnm) )
1535 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1541 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1544 "\nWarning: assuming the time unit in %s is %s\n",
1545 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1547 if (trx_out_fn && !bReadTraj)
1549 fprintf(stderr, "\nWarning: "
1550 "cannot write cluster structures without reading trajectory\n"
1551 " ignoring option -cl %s\n", trx_out_fn);
1555 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1561 gmx_fatal(FARGS, "Invalid method");
1564 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1565 method == m_gromos );
1568 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1570 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1571 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1573 /* check input and write parameters to log file */
1574 bUseRmsdCut = FALSE;
1575 if (method == m_jarvis_patrick)
1577 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1578 if ((M < 0) || (M == 1))
1580 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1584 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1591 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1595 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1600 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1603 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1605 else /* method != m_jarvis */
1607 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1609 if (bUseRmsdCut && method != m_jarvis_patrick)
1611 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1613 if (method == m_monte_carlo)
1615 fprintf(log, "Using %d iterations\n", niter);
1620 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1626 /* don't read mass-database as masses (and top) are not used */
1627 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1631 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1634 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1635 bReadMat ? "" : " and RMSD calculation");
1636 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1637 1, &ifsize, &fitidx, &grpname);
1640 fprintf(stderr, "\nSelect group for output:\n");
1641 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1642 1, &iosize, &outidx, &grpname);
1643 /* merge and convert both index groups: */
1644 /* first copy outidx to index. let outidx refer to elements in index */
1645 snew(index, iosize);
1647 for (i = 0; i < iosize; i++)
1649 index[i] = outidx[i];
1652 /* now lookup elements from fitidx in index, add them if necessary
1653 and also let fitidx refer to elements in index */
1654 for (i = 0; i < ifsize; i++)
1657 while (j < isize && index[j] != fitidx[i])
1663 /* slow this way, but doesn't matter much */
1665 srenew(index, isize);
1667 index[j] = fitidx[i];
1671 else /* !trx_out_fn */
1675 for (i = 0; i < ifsize; i++)
1677 index[i] = fitidx[i];
1685 /* Loop over first coordinate file */
1686 fn = opt2fn("-f", NFILE, fnm);
1688 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1689 output_env_conv_times(oenv, nf, time);
1690 if (!bRMSdist || bAnalyze)
1692 /* Center all frames on zero */
1694 for (i = 0; i < ifsize; i++)
1696 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1700 for (i = 0; i < nf; i++)
1702 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1708 gmx_rmpbc_done(gpbc);
1714 fprintf(stderr, "Reading rms distance matrix ");
1715 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1716 fprintf(stderr, "\n");
1717 if (readmat[0].nx != readmat[0].ny)
1719 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1720 readmat[0].nx, readmat[0].ny);
1722 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1724 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1725 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1730 time = readmat[0].axis_x;
1731 time_invfac = output_env_get_time_invfactor(oenv);
1732 for (i = 0; i < nf; i++)
1734 time[i] *= time_invfac;
1737 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1738 convert_mat(&(readmat[0]), rms);
1740 nlevels = readmat[0].nmap;
1742 else /* !bReadMat */
1744 rms = init_mat(nf, method == m_diagonalize);
1745 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1748 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1749 /* Initialize work array */
1751 for (i1 = 0; i1 < nf; i1++)
1753 for (i2 = i1+1; i2 < nf; i2++)
1755 for (i = 0; i < isize; i++)
1757 copy_rvec(xx[i1][i], x1[i]);
1761 do_fit(isize, mass, xx[i2], x1);
1763 rmsd = rmsdev(isize, mass, xx[i2], x1);
1764 set_mat_entry(rms, i1, i2, rmsd);
1766 nrms -= (gmx_int64_t) (nf-i1-1);
1767 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1773 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1775 /* Initiate work arrays */
1778 for (i = 0; (i < isize); i++)
1783 for (i1 = 0; i1 < nf; i1++)
1785 calc_dist(isize, xx[i1], d1);
1786 for (i2 = i1+1; (i2 < nf); i2++)
1788 calc_dist(isize, xx[i2], d2);
1789 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1792 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1794 /* Clean up work arrays */
1795 for (i = 0; (i < isize); i++)
1803 fprintf(stderr, "\n\n");
1805 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1806 rms->minrms, rms->maxrms);
1807 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1808 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1809 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1810 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1812 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1813 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1815 if (bAnalyze && (rmsmin < rms->minrms) )
1817 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1818 rmsmin, rms->minrms);
1820 if (bAnalyze && (rmsmin > rmsdcut) )
1822 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1826 /* Plot the rmsd distribution */
1827 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1831 for (i1 = 0; (i1 < nf); i1++)
1833 for (i2 = 0; (i2 < nf); i2++)
1835 if (rms->mat[i1][i2] < rmsdcut)
1837 rms->mat[i1][i2] = 0;
1841 rms->mat[i1][i2] = 1;
1851 /* Now sort the matrix and write it out again */
1852 gather(rms, rmsdcut, &clust);
1855 /* Do a diagonalization */
1856 snew(eigenvalues, nf);
1857 snew(eigenvectors, nf*nf);
1858 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1859 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1860 sfree(eigenvectors);
1862 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1863 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1864 for (i = 0; (i < nf); i++)
1866 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1871 orig = init_mat(rms->nn, FALSE);
1873 copy_t_mat(orig, rms);
1874 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1875 opt2fn_null("-conv", NFILE, fnm), oenv);
1877 case m_jarvis_patrick:
1878 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1881 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1884 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1887 if (method == m_monte_carlo || method == m_diagonalize)
1889 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1897 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1901 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1903 init_t_atoms(&useatoms, isize, FALSE);
1904 snew(usextps, isize);
1905 useatoms.resinfo = top.atoms.resinfo;
1906 for (i = 0; i < isize; i++)
1908 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1909 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1910 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1911 copy_rvec(xtps[index[i]], usextps[i]);
1913 useatoms.nr = isize;
1914 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1915 ifsize, fitidx, iosize, outidx,
1916 bReadTraj ? trx_out_fn : NULL,
1917 opt2fn_null("-sz", NFILE, fnm),
1918 opt2fn_null("-tr", NFILE, fnm),
1919 opt2fn_null("-ntr", NFILE, fnm),
1920 opt2fn_null("-clid", NFILE, fnm),
1921 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1922 rlo_bot, rhi_bot, oenv);
1926 if (bBinary && !bAnalyze)
1928 /* Make the clustering visible */
1929 for (i2 = 0; (i2 < nf); i2++)
1931 for (i1 = i2+1; (i1 < nf); i1++)
1933 if (rms->mat[i1][i2])
1935 rms->mat[i1][i2] = rms->maxrms;
1941 fp = opt2FILE("-o", NFILE, fnm, "w");
1942 fprintf(stderr, "Writing rms distance/clustering matrix ");
1945 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1946 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1947 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1951 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1952 sprintf(title, "RMS%sDeviation / Cluster Index",
1953 bRMSdist ? " Distance " : " ");
1956 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1957 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1958 rlo_top, rhi_top, 0.0, (real) ncluster,
1959 &ncluster, TRUE, rlo_bot, rhi_bot);
1963 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1964 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1965 rlo_top, rhi_top, &nlevels);
1968 fprintf(stderr, "\n");
1972 fp = opt2FILE("-om", NFILE, fnm, "w");
1973 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1974 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1975 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1976 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1977 rlo_top, rhi_top, &nlevels);
1982 /* now show what we've done */
1983 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1984 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1985 if (method == m_diagonalize)
1987 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1989 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1992 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1993 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1994 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1996 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);