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.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/tpxio.h"
49 #include "gromacs/fileio/trxio.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"
68 /* print to two file pointers at once (i.e. stderr and log) */
70 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
72 fprintf(fp1, "%s", buf);
73 fprintf(fp2, "%s", buf);
76 /* just print a prepared buffer to fp1 and fp2 */
78 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
80 lo_ffprintf(fp1, fp2, buf);
83 /* prepare buffer with one argument, then print to fp1 and fp2 */
85 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
87 sprintf(buf, fmt, arg);
88 lo_ffprintf(fp1, fp2, buf);
91 /* prepare buffer with one argument, then print to fp1 and fp2 */
93 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
95 sprintf(buf, fmt, arg);
96 lo_ffprintf(fp1, fp2, buf);
99 /* prepare buffer with one argument, then print to fp1 and fp2 */
101 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
103 sprintf(buf, fmt, arg);
104 lo_ffprintf(fp1, fp2, buf);
107 /* prepare buffer with two arguments, then print to fp1 and fp2 */
109 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
111 sprintf(buf, fmt, arg1, arg2);
112 lo_ffprintf(fp1, fp2, buf);
115 /* prepare buffer with two arguments, then print to fp1 and fp2 */
117 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
119 sprintf(buf, fmt, arg1, arg2);
120 lo_ffprintf(fp1, fp2, buf);
123 /* prepare buffer with two arguments, then print to fp1 and fp2 */
125 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
127 sprintf(buf, fmt, arg1, arg2);
128 lo_ffprintf(fp1, fp2, buf);
141 void pr_energy(FILE *fp, real e)
143 fprintf(fp, "Energy: %8.4f\n", e);
146 void cp_index(int nn, int from[], int to[])
150 for (i = 0; (i < nn); i++)
156 void mc_optimize(FILE *log, t_mat *m, real *time,
157 int maxiter, int nrandom,
159 const char *conv, output_env_t oenv)
162 real ecur, enext, emin, prob, enorm;
163 int i, j, iswap, jswap, nn, nuphill = 0;
169 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
172 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
173 printf("by reordering the frames to minimize the path between the two structures\n");
174 printf("that have the largest pairwise RMSD.\n");
177 enorm = m->mat[0][0];
178 for (i = 0; (i < m->n1); i++)
180 for (j = 0; (j < m->nn); j++)
182 if (m->mat[i][j] > enorm)
184 enorm = m->mat[i][j];
190 if ((iswap == -1) || (jswap == -1))
192 fprintf(stderr, "Matrix contains identical values in all fields\n");
195 swap_rows(m, 0, iswap);
196 swap_rows(m, m->n1-1, jswap);
197 emin = ecur = mat_energy(m);
198 printf("Largest distance %g between %d and %d. Energy: %g.\n",
199 enorm, iswap, jswap, emin);
201 rng = gmx_rng_init(seed);
204 /* Initiate and store global minimum */
205 minimum = init_mat(nn, m->b1D);
207 copy_t_mat(minimum, m);
211 fp = xvgropen(conv, "Convergence of the MC optimization",
212 "Energy", "Step", oenv);
214 for (i = 0; (i < maxiter); i++)
216 /* Generate new swapping candidates */
219 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
220 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
222 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
224 /* Apply swap and compute energy */
225 swap_rows(m, iswap, jswap);
226 enext = mat_energy(m);
228 /* Compute probability */
230 if ((enext < ecur) || (i < nrandom))
235 /* Store global minimum */
236 copy_t_mat(minimum, m);
242 /* Try Monte Carlo step */
243 prob = exp(-(enext-ecur)/(enorm*kT));
246 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
253 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
254 i, iswap, jswap, enext, prob);
257 fprintf(fp, "%6d %10g\n", i, enext);
263 swap_rows(m, jswap, iswap);
266 fprintf(log, "%d uphill steps were taken during optimization\n",
269 /* Now swap the matrix to get it into global minimum mode */
270 copy_t_mat(m, minimum);
272 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
273 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
274 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
275 for (i = 0; (i < m->nn); i++)
277 fprintf(log, "%10g %5d %10g\n",
280 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
289 static void calc_dist(int nind, rvec x[], real **d)
295 for (i = 0; (i < nind-1); i++)
298 for (j = i+1; (j < nind); j++)
300 /* Should use pbc_dx when analysing multiple molecueles,
301 * but the box is not stored for every frame.
303 rvec_sub(xi, x[j], dx);
309 static real rms_dist(int isize, real **d, real **d_r)
315 for (i = 0; (i < isize-1); i++)
317 for (j = i+1; (j < isize); j++)
319 r = d[i][j]-d_r[i][j];
323 r2 /= (isize*(isize-1))/2;
328 static int rms_dist_comp(const void *a, const void *b)
335 if (da->dist - db->dist < 0)
339 else if (da->dist - db->dist > 0)
346 static int clust_id_comp(const void *a, const void *b)
353 return da->clust - db->clust;
356 static int nrnb_comp(const void *a, const void *b)
363 /* return the b-a, we want highest first */
364 return db->nr - da->nr;
367 void gather(t_mat *m, real cutoff, t_clusters *clust)
371 int i, j, k, nn, cid, n1, diff;
374 /* First we sort the entries in the RMSD matrix */
378 for (i = k = 0; (i < n1); i++)
380 for (j = i+1; (j < n1); j++, k++)
384 d[k].dist = m->mat[i][j];
389 gmx_incons("gather algortihm");
391 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
393 /* Now we make a cluster index for all of the conformations */
396 /* Now we check the closest structures, and equalize their cluster numbers */
397 fprintf(stderr, "Linking structures ");
400 fprintf(stderr, "*");
402 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
404 diff = c[d[k].j].clust - c[d[k].i].clust;
410 c[d[k].j].clust = c[d[k].i].clust;
414 c[d[k].i].clust = c[d[k].j].clust;
420 fprintf(stderr, "\nSorting and renumbering clusters\n");
421 /* Sort on cluster number */
422 qsort(c, n1, sizeof(c[0]), clust_id_comp);
424 /* Renumber clusters */
426 for (k = 1; k < n1; k++)
428 if (c[k].clust != c[k-1].clust)
441 for (k = 0; (k < n1); k++)
443 fprintf(debug, "Cluster index for conformation %d: %d\n",
444 c[k].conf, c[k].clust);
448 for (k = 0; k < n1; k++)
450 clust->cl[c[k].conf] = c[k].clust;
457 gmx_bool jp_same(int **nnb, int i, int j, int P)
463 for (k = 0; nnb[i][k] >= 0; k++)
465 bIn = bIn || (nnb[i][k] == j);
473 for (k = 0; nnb[j][k] >= 0; k++)
475 bIn = bIn || (nnb[j][k] == i);
483 for (ii = 0; nnb[i][ii] >= 0; ii++)
485 for (jj = 0; nnb[j][jj] >= 0; jj++)
487 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
497 static void jarvis_patrick(int n1, real **mat, int M, int P,
498 real rmsdcut, t_clusters *clust)
503 int i, j, k, cid, diff, max;
512 /* First we sort the entries in the RMSD matrix row by row.
513 * This gives us the nearest neighbor list.
517 for (i = 0; (i < n1); i++)
519 for (j = 0; (j < n1); j++)
522 row[j].dist = mat[i][j];
524 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
527 /* Put the M nearest neighbors in the list */
529 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
533 nnb[i][k] = row[j].j;
541 /* Put all neighbors nearer than rmsdcut in the list */
544 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
553 nnb[i][k] = row[j].j;
559 srenew(nnb[i], max+1);
567 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
568 for (i = 0; (i < n1); i++)
570 fprintf(debug, "i:%5d nbs:", i);
571 for (j = 0; nnb[i][j] >= 0; j++)
573 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
575 fprintf(debug, "\n");
580 fprintf(stderr, "Linking structures ");
581 /* Use mcpy for temporary storage of booleans */
582 mcpy = mk_matrix(n1, n1, FALSE);
583 for (i = 0; i < n1; i++)
585 for (j = i+1; j < n1; j++)
587 mcpy[i][j] = jp_same(nnb, i, j, P);
592 fprintf(stderr, "*");
594 for (i = 0; i < n1; i++)
596 for (j = i+1; j < n1; j++)
600 diff = c[j].clust - c[i].clust;
606 c[j].clust = c[i].clust;
610 c[i].clust = c[j].clust;
619 fprintf(stderr, "\nSorting and renumbering clusters\n");
620 /* Sort on cluster number */
621 qsort(c, n1, sizeof(c[0]), clust_id_comp);
623 /* Renumber clusters */
625 for (k = 1; k < n1; k++)
627 if (c[k].clust != c[k-1].clust)
639 for (k = 0; k < n1; k++)
641 clust->cl[c[k].conf] = c[k].clust;
645 for (k = 0; (k < n1); k++)
647 fprintf(debug, "Cluster index for conformation %d: %d\n",
648 c[k].conf, c[k].clust);
652 /* Again, I don't see the point in this... (AF) */
653 /* for(i=0; (i<n1); i++) { */
654 /* for(j=0; (j<n1); j++) */
655 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
657 /* for(i=0; (i<n1); i++) { */
658 /* for(j=0; (j<n1); j++) */
659 /* mat[i][j] = mcpy[i][j]; */
661 done_matrix(n1, &mcpy);
664 for (i = 0; (i < n1); i++)
671 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
675 /* dump neighbor list */
676 fprintf(fp, "%s", title);
677 for (i = 0; (i < n1); i++)
679 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
680 for (j = 0; j < nnb[i].nr; j++)
682 fprintf(fp, "%5d", nnb[i].nb[j]);
688 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
692 int i, j, k, j1, max;
694 /* Put all neighbors nearer than rmsdcut in the list */
695 fprintf(stderr, "Making list of neighbors within cutoff ");
698 for (i = 0; (i < n1); i++)
702 /* put all neighbors within cut-off in list */
703 for (j = 0; j < n1; j++)
705 if (mat[i][j] < rmsdcut)
710 srenew(nnb[i].nb, max);
716 /* store nr of neighbors, we'll need that */
718 if (i%(1+n1/100) == 0)
720 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
723 fprintf(stderr, "%3d%%\n", 100);
726 /* sort neighbor list on number of neighbors, largest first */
727 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
731 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
734 /* turn first structure with all its neighbors (largest) into cluster
735 remove them from pool of structures and repeat for all remaining */
736 fprintf(stderr, "Finding clusters %4d", 0);
737 /* cluster id's start at 1: */
741 /* set cluster id (k) for first item in neighborlist */
742 for (j = 0; j < nnb[0].nr; j++)
744 clust->cl[nnb[0].nb[j]] = k;
750 /* adjust number of neighbors for others, taking removals into account: */
751 for (i = 1; i < n1 && nnb[i].nr; i++)
754 for (j = 0; j < nnb[i].nr; j++)
756 /* if this neighbor wasn't removed */
757 if (clust->cl[nnb[i].nb[j]] == 0)
759 /* shift the rest (j1<=j) */
760 nnb[i].nb[j1] = nnb[i].nb[j];
765 /* now j1 is the new number of neighbors */
768 /* sort again on nnb[].nr, because we have new # neighbors: */
769 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
770 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
772 fprintf(stderr, "\b\b\b\b%4d", k);
776 fprintf(stderr, "\n");
780 fprintf(debug, "Clusters (%d):\n", k);
781 for (i = 0; i < n1; i++)
783 fprintf(debug, " %3d", clust->cl[i]);
785 fprintf(debug, "\n");
791 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
792 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
797 int i, i0, j, max_nf;
805 natom = read_first_x(oenv, &status, fn, &t, &x, box);
812 gmx_rmpbc(gpbc, natom, box, x);
818 srenew(*time, max_nf);
823 /* Store only the interesting atoms */
824 for (j = 0; (j < isize); j++)
826 copy_rvec(x[index[j]], xx[i0][j]);
833 while (read_next_x(oenv, status, &t, x, box));
834 fprintf(stderr, "Allocated %lu bytes for frames\n",
835 (unsigned long) (max_nf*isize*sizeof(**xx)));
836 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
843 static int plot_clusters(int nf, real **mat, t_clusters *clust,
846 int i, j, ncluster, ci;
847 int *cl_id, *nstruct, *strind;
852 for (i = 0; i < nf; i++)
855 cl_id[i] = clust->cl[i];
859 for (i = 0; i < nf; i++)
861 if (nstruct[i] >= minstruct)
864 for (j = 0; (j < nf); j++)
868 strind[j] = ncluster;
874 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
875 ncluster, minstruct);
877 for (i = 0; (i < nf); i++)
880 for (j = 0; j < i; j++)
882 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
884 /* color different clusters with different colors, as long as
885 we don't run out of colors */
886 mat[i][j] = strind[i];
901 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
905 for (i = 0; i < nf; i++)
907 for (j = 0; j < i; j++)
909 if (clust->cl[i] == clust->cl[j])
921 static char *parse_filename(const char *fn, int maxnr)
930 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
932 /* number of digits needed in numbering */
933 i = (int)(log(maxnr)/log(10)) + 1;
934 /* split fn and ext */
935 ext = strrchr(fn, '.');
938 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
941 /* insert e.g. '%03d' between fn and ext */
942 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
943 snew(fnout, strlen(buf)+1);
949 static void ana_trans(t_clusters *clust, int nf,
950 const char *transfn, const char *ntransfn, FILE *log,
951 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
956 int i, ntranst, maxtrans;
959 snew(ntrans, clust->ncl);
960 snew(trans, clust->ncl);
961 snew(axis, clust->ncl);
962 for (i = 0; i < clust->ncl; i++)
965 snew(trans[i], clust->ncl);
969 for (i = 1; i < nf; i++)
971 if (clust->cl[i] != clust->cl[i-1])
974 ntrans[clust->cl[i-1]-1]++;
975 ntrans[clust->cl[i]-1]++;
976 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
977 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
980 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
981 "max %d between two specific clusters\n", ntranst, maxtrans);
984 fp = gmx_ffopen(transfn, "w");
985 i = min(maxtrans+1, 80);
986 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
987 "from cluster", "to cluster",
988 clust->ncl, clust->ncl, axis, axis, trans,
989 0, maxtrans, rlo, rhi, &i);
994 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
996 for (i = 0; i < clust->ncl; i++)
998 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1003 for (i = 0; i < clust->ncl; i++)
1011 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1012 int natom, t_atoms *atoms, rvec *xtps,
1013 real *mass, rvec **xx, real *time,
1014 int ifsize, atom_id *fitidx,
1015 int iosize, atom_id *outidx,
1016 const char *trxfn, const char *sizefn,
1017 const char *transfn, const char *ntransfn,
1018 const char *clustidfn, gmx_bool bAverage,
1019 int write_ncl, int write_nst, real rmsmin,
1020 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1021 const output_env_t oenv)
1024 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1025 t_trxstatus *trxout = NULL;
1026 t_trxstatus *trxsout = NULL;
1027 int i, i1, cl, nstr, *structure, first = 0, midstr;
1028 gmx_bool *bWrite = NULL;
1029 real r, clrmsd, midrmsd;
1035 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1039 /* do we write all structures? */
1042 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1045 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1046 bAverage ? "average" : "middle", trxfn);
1049 /* find out what we want to tell the user:
1050 Writing [all structures|structures with rmsd > %g] for
1051 {all|first %d} clusters {with more than %d structures} to %s */
1054 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1058 sprintf(buf1, "all structures");
1060 buf2[0] = buf3[0] = '\0';
1061 if (write_ncl >= clust->ncl)
1065 sprintf(buf2, "all ");
1070 sprintf(buf2, "the first %d ", write_ncl);
1074 sprintf(buf3, " with more than %d structures", write_nst);
1076 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1077 ffprintf(stderr, log, buf);
1080 /* Prepare a reference structure for the orientation of the clusters */
1083 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1085 trxout = open_trx(trxfn, "w");
1086 /* Calculate the average structure in each cluster, *
1087 * all structures are fitted to the first struture of the cluster */
1091 if (transfn || ntransfn)
1093 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1098 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1099 fprintf(fp, "@ s0 symbol 2\n");
1100 fprintf(fp, "@ s0 symbol size 0.2\n");
1101 fprintf(fp, "@ s0 linestyle 0\n");
1102 for (i = 0; i < nf; i++)
1104 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1110 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1111 fprintf(fp, "@g%d type %s\n", 0, "bar");
1113 snew(structure, nf);
1114 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1115 "cl.", "#st", "rmsd", "middle", "rmsd");
1116 for (cl = 1; cl <= clust->ncl; cl++)
1118 /* prepare structures (fit, middle, average) */
1121 for (i = 0; i < natom; i++)
1127 for (i1 = 0; i1 < nf; i1++)
1129 if (clust->cl[i1] == cl)
1131 structure[nstr] = i1;
1133 if (trxfn && (bAverage || write_ncl) )
1137 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1145 do_fit(natom, mass, xx[first], xx[i1]);
1149 for (i = 0; i < natom; i++)
1151 rvec_inc(xav[i], xx[i1][i]);
1159 fprintf(fp, "%8d %8d\n", cl, nstr);
1164 for (i1 = 0; i1 < nstr; i1++)
1169 for (i = 0; i < nstr; i++)
1173 r += rmsd[structure[i]][structure[i1]];
1177 r += rmsd[structure[i1]][structure[i]];
1184 midstr = structure[i1];
1191 /* dump cluster info to logfile */
1194 sprintf(buf1, "%6.3f", clrmsd);
1199 sprintf(buf2, "%5.3f", midrmsd);
1207 sprintf(buf1, "%5s", "");
1208 sprintf(buf2, "%5s", "");
1210 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1211 for (i = 0; i < nstr; i++)
1213 if ((i % 7 == 0) && i)
1215 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1222 fprintf(log, "%s %6g", buf, time[i1]);
1226 /* write structures to trajectory file(s) */
1231 for (i = 0; i < nstr; i++)
1236 if (cl < write_ncl+1 && nstr > write_nst)
1238 /* Dump all structures for this cluster */
1239 /* generate numbered filename (there is a %d in trxfn!) */
1240 sprintf(buf, trxsfn, cl);
1241 trxsout = open_trx(buf, "w");
1242 for (i = 0; i < nstr; i++)
1247 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1251 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1257 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1258 xx[structure[i]], NULL, NULL);
1263 /* Dump the average structure for this cluster */
1266 for (i = 0; i < natom; i++)
1268 svmul(1.0/nstr, xav[i], xav[i]);
1273 for (i = 0; i < natom; i++)
1275 copy_rvec(xx[midstr][i], xav[i]);
1279 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1284 do_fit(natom, mass, xtps, xav);
1287 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1307 static void convert_mat(t_matrix *mat, t_mat *rms)
1312 matrix2real(mat, rms->mat);
1313 /* free input xpm matrix data */
1314 for (i = 0; i < mat->nx; i++)
1316 sfree(mat->matrix[i]);
1320 for (i = 0; i < mat->nx; i++)
1322 for (j = i; j < mat->nx; j++)
1324 rms->sumrms += rms->mat[i][j];
1325 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1328 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1335 int gmx_cluster(int argc, char *argv[])
1337 const char *desc[] = {
1338 "[THISMODULE] can cluster structures using several different methods.",
1339 "Distances between structures can be determined from a trajectory",
1340 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1341 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1342 "can be used to define the distance between structures.[PAR]",
1344 "single linkage: add a structure to a cluster when its distance to any",
1345 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1347 "Jarvis Patrick: add a structure to a cluster when this structure",
1348 "and a structure in the cluster have each other as neighbors and",
1349 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1350 "of a structure are the M closest structures or all structures within",
1351 "[TT]cutoff[tt].[PAR]",
1353 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1354 "the order of the frames is using the smallest possible increments.",
1355 "With this it is possible to make a smooth animation going from one",
1356 "structure to another with the largest possible (e.g.) RMSD between",
1357 "them, however the intermediate steps should be as small as possible.",
1358 "Applications could be to visualize a potential of mean force",
1359 "ensemble of simulations or a pulling simulation. Obviously the user",
1360 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1361 "The final result can be inspect visually by looking at the matrix",
1362 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1364 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1366 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1367 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1368 "Count number of neighbors using cut-off, take structure with",
1369 "largest number of neighbors with all its neighbors as cluster",
1370 "and eliminate it from the pool of clusters. Repeat for remaining",
1371 "structures in pool.[PAR]",
1373 "When the clustering algorithm assigns each structure to exactly one",
1374 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1375 "file is supplied, the structure with",
1376 "the smallest average distance to the others or the average structure",
1377 "or all structures for each cluster will be written to a trajectory",
1378 "file. When writing all structures, separate numbered files are made",
1379 "for each cluster.[PAR]",
1381 "Two output files are always written:[BR]",
1382 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1383 "and a graphical depiction of the clusters in the lower right half",
1384 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1385 "when two structures are in the same cluster.",
1386 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1388 "[TT]-g[tt] writes information on the options used and a detailed list",
1389 "of all clusters and their members.[PAR]",
1391 "Additionally, a number of optional output files can be written:[BR]",
1392 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1393 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1394 "diagonalization.[BR]",
1395 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1396 "[TT]-tr[tt] writes a matrix of the number transitions between",
1397 "cluster pairs.[BR]",
1398 "[TT]-ntr[tt] writes the total number of transitions to or from",
1399 "each cluster.[BR]",
1400 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1401 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1402 "structure of each cluster or writes numbered files with cluster members",
1403 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1404 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1405 "structure with the smallest average RMSD from all other structures",
1406 "of the cluster.[BR]",
1410 int nf, i, i1, i2, j;
1411 gmx_int64_t nrms = 0;
1414 rvec *xtps, *usextps, *x1, **xx = NULL;
1415 const char *fn, *trx_out_fn;
1417 t_mat *rms, *orig = NULL;
1422 t_matrix *readmat = NULL;
1425 int isize = 0, ifsize = 0, iosize = 0;
1426 atom_id *index = NULL, *fitidx, *outidx;
1428 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1429 char buf[STRLEN], buf1[80], title[STRLEN];
1430 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1432 int method, ncluster = 0;
1433 static const char *methodname[] = {
1434 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1435 "diagonalization", "gromos", NULL
1438 m_null, m_linkage, m_jarvis_patrick,
1439 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1441 /* Set colors for plotting: white = zero RMS, black = maximum */
1442 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1443 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1444 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1445 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1446 static int nlevels = 40, skip = 1;
1447 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1448 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1449 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1450 static real kT = 1e-3;
1451 static int M = 10, P = 3;
1453 gmx_rmpbc_t gpbc = NULL;
1456 { "-dista", FALSE, etBOOL, {&bRMSdist},
1457 "Use RMSD of distances instead of RMS deviation" },
1458 { "-nlevels", FALSE, etINT, {&nlevels},
1459 "Discretize RMSD matrix in this number of levels" },
1460 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1461 "RMSD cut-off (nm) for two structures to be neighbor" },
1462 { "-fit", FALSE, etBOOL, {&bFit},
1463 "Use least squares fitting before RMSD calculation" },
1464 { "-max", FALSE, etREAL, {&scalemax},
1465 "Maximum level in RMSD matrix" },
1466 { "-skip", FALSE, etINT, {&skip},
1467 "Only analyze every nr-th frame" },
1468 { "-av", FALSE, etBOOL, {&bAverage},
1469 "Write average iso middle structure for each cluster" },
1470 { "-wcl", FALSE, etINT, {&write_ncl},
1471 "Write the structures for this number of clusters to numbered files" },
1472 { "-nst", FALSE, etINT, {&write_nst},
1473 "Only write all structures if more than this number of structures per cluster" },
1474 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1475 "minimum rms difference with rest of cluster for writing structures" },
1476 { "-method", FALSE, etENUM, {methodname},
1477 "Method for cluster determination" },
1478 { "-minstruct", FALSE, etINT, {&minstruct},
1479 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1480 { "-binary", FALSE, etBOOL, {&bBinary},
1481 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1482 "is given by [TT]-cutoff[tt]" },
1483 { "-M", FALSE, etINT, {&M},
1484 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1485 "0 is use cutoff" },
1486 { "-P", FALSE, etINT, {&P},
1487 "Number of identical nearest neighbors required to form a cluster" },
1488 { "-seed", FALSE, etINT, {&seed},
1489 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1490 { "-niter", FALSE, etINT, {&niter},
1491 "Number of iterations for MC" },
1492 { "-nrandom", FALSE, etINT, {&nrandom},
1493 "The first iterations for MC may be done complete random, to shuffle the frames" },
1494 { "-kT", FALSE, etREAL, {&kT},
1495 "Boltzmann weighting factor for Monte Carlo optimization "
1496 "(zero turns off uphill steps)" },
1497 { "-pbc", FALSE, etBOOL,
1498 { &bPBC }, "PBC check" }
1501 { efTRX, "-f", NULL, ffOPTRD },
1502 { efTPS, "-s", NULL, ffOPTRD },
1503 { efNDX, NULL, NULL, ffOPTRD },
1504 { efXPM, "-dm", "rmsd", ffOPTRD },
1505 { efXPM, "-om", "rmsd-raw", ffWRITE },
1506 { efXPM, "-o", "rmsd-clust", ffWRITE },
1507 { efLOG, "-g", "cluster", ffWRITE },
1508 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1509 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1510 { efXVG, "-conv", "mc-conv", ffOPTWR },
1511 { efXVG, "-sz", "clust-size", ffOPTWR},
1512 { efXPM, "-tr", "clust-trans", ffOPTWR},
1513 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1514 { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1515 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1517 #define NFILE asize(fnm)
1519 if (!parse_common_args(&argc, argv,
1520 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1521 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1528 bReadMat = opt2bSet("-dm", NFILE, fnm);
1529 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1530 if (opt2parg_bSet("-av", asize(pa), pa) ||
1531 opt2parg_bSet("-wcl", asize(pa), pa) ||
1532 opt2parg_bSet("-nst", asize(pa), pa) ||
1533 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1534 opt2bSet("-cl", NFILE, fnm) )
1536 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1542 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1545 "\nWarning: assuming the time unit in %s is %s\n",
1546 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1548 if (trx_out_fn && !bReadTraj)
1550 fprintf(stderr, "\nWarning: "
1551 "cannot write cluster structures without reading trajectory\n"
1552 " ignoring option -cl %s\n", trx_out_fn);
1556 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1562 gmx_fatal(FARGS, "Invalid method");
1565 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1566 method == m_gromos );
1569 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1571 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1572 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1574 /* check input and write parameters to log file */
1575 bUseRmsdCut = FALSE;
1576 if (method == m_jarvis_patrick)
1578 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1579 if ((M < 0) || (M == 1))
1581 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1585 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1592 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1596 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1601 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1604 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1606 else /* method != m_jarvis */
1608 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1610 if (bUseRmsdCut && method != m_jarvis_patrick)
1612 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1614 if (method == m_monte_carlo)
1616 fprintf(log, "Using %d iterations\n", niter);
1621 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1627 /* don't read mass-database as masses (and top) are not used */
1628 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1632 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1635 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1636 bReadMat ? "" : " and RMSD calculation");
1637 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1638 1, &ifsize, &fitidx, &grpname);
1641 fprintf(stderr, "\nSelect group for output:\n");
1642 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1643 1, &iosize, &outidx, &grpname);
1644 /* merge and convert both index groups: */
1645 /* first copy outidx to index. let outidx refer to elements in index */
1646 snew(index, iosize);
1648 for (i = 0; i < iosize; i++)
1650 index[i] = outidx[i];
1653 /* now lookup elements from fitidx in index, add them if necessary
1654 and also let fitidx refer to elements in index */
1655 for (i = 0; i < ifsize; i++)
1658 while (j < isize && index[j] != fitidx[i])
1664 /* slow this way, but doesn't matter much */
1666 srenew(index, isize);
1668 index[j] = fitidx[i];
1672 else /* !trx_out_fn */
1676 for (i = 0; i < ifsize; i++)
1678 index[i] = fitidx[i];
1686 /* Loop over first coordinate file */
1687 fn = opt2fn("-f", NFILE, fnm);
1689 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1690 output_env_conv_times(oenv, nf, time);
1691 if (!bRMSdist || bAnalyze)
1693 /* Center all frames on zero */
1695 for (i = 0; i < ifsize; i++)
1697 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1701 for (i = 0; i < nf; i++)
1703 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1709 gmx_rmpbc_done(gpbc);
1715 fprintf(stderr, "Reading rms distance matrix ");
1716 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1717 fprintf(stderr, "\n");
1718 if (readmat[0].nx != readmat[0].ny)
1720 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1721 readmat[0].nx, readmat[0].ny);
1723 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1725 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1726 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1731 time = readmat[0].axis_x;
1732 time_invfac = output_env_get_time_invfactor(oenv);
1733 for (i = 0; i < nf; i++)
1735 time[i] *= time_invfac;
1738 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1739 convert_mat(&(readmat[0]), rms);
1741 nlevels = readmat[0].nmap;
1743 else /* !bReadMat */
1745 rms = init_mat(nf, method == m_diagonalize);
1746 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1749 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1750 /* Initialize work array */
1752 for (i1 = 0; i1 < nf; i1++)
1754 for (i2 = i1+1; i2 < nf; i2++)
1756 for (i = 0; i < isize; i++)
1758 copy_rvec(xx[i1][i], x1[i]);
1762 do_fit(isize, mass, xx[i2], x1);
1764 rmsd = rmsdev(isize, mass, xx[i2], x1);
1765 set_mat_entry(rms, i1, i2, rmsd);
1767 nrms -= (gmx_int64_t) (nf-i1-1);
1768 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1774 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1776 /* Initiate work arrays */
1779 for (i = 0; (i < isize); i++)
1784 for (i1 = 0; i1 < nf; i1++)
1786 calc_dist(isize, xx[i1], d1);
1787 for (i2 = i1+1; (i2 < nf); i2++)
1789 calc_dist(isize, xx[i2], d2);
1790 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1793 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1795 /* Clean up work arrays */
1796 for (i = 0; (i < isize); i++)
1804 fprintf(stderr, "\n\n");
1806 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1807 rms->minrms, rms->maxrms);
1808 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1809 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1810 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1811 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1813 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1814 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1816 if (bAnalyze && (rmsmin < rms->minrms) )
1818 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1819 rmsmin, rms->minrms);
1821 if (bAnalyze && (rmsmin > rmsdcut) )
1823 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1827 /* Plot the rmsd distribution */
1828 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1832 for (i1 = 0; (i1 < nf); i1++)
1834 for (i2 = 0; (i2 < nf); i2++)
1836 if (rms->mat[i1][i2] < rmsdcut)
1838 rms->mat[i1][i2] = 0;
1842 rms->mat[i1][i2] = 1;
1852 /* Now sort the matrix and write it out again */
1853 gather(rms, rmsdcut, &clust);
1856 /* Do a diagonalization */
1857 snew(eigenvalues, nf);
1858 snew(eigenvectors, nf*nf);
1859 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1860 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1861 sfree(eigenvectors);
1863 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1864 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1865 for (i = 0; (i < nf); i++)
1867 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1872 orig = init_mat(rms->nn, FALSE);
1874 copy_t_mat(orig, rms);
1875 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1876 opt2fn_null("-conv", NFILE, fnm), oenv);
1878 case m_jarvis_patrick:
1879 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1882 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1885 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1888 if (method == m_monte_carlo || method == m_diagonalize)
1890 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1898 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1902 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1904 init_t_atoms(&useatoms, isize, FALSE);
1905 snew(usextps, isize);
1906 useatoms.resinfo = top.atoms.resinfo;
1907 for (i = 0; i < isize; i++)
1909 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1910 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1911 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1912 copy_rvec(xtps[index[i]], usextps[i]);
1914 useatoms.nr = isize;
1915 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1916 ifsize, fitidx, iosize, outidx,
1917 bReadTraj ? trx_out_fn : NULL,
1918 opt2fn_null("-sz", NFILE, fnm),
1919 opt2fn_null("-tr", NFILE, fnm),
1920 opt2fn_null("-ntr", NFILE, fnm),
1921 opt2fn_null("-clid", NFILE, fnm),
1922 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1923 rlo_bot, rhi_bot, oenv);
1927 if (bBinary && !bAnalyze)
1929 /* Make the clustering visible */
1930 for (i2 = 0; (i2 < nf); i2++)
1932 for (i1 = i2+1; (i1 < nf); i1++)
1934 if (rms->mat[i1][i2])
1936 rms->mat[i1][i2] = rms->maxrms;
1942 fp = opt2FILE("-o", NFILE, fnm, "w");
1943 fprintf(stderr, "Writing rms distance/clustering matrix ");
1946 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1947 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1948 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1952 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1953 sprintf(title, "RMS%sDeviation / Cluster Index",
1954 bRMSdist ? " Distance " : " ");
1957 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1958 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1959 rlo_top, rhi_top, 0.0, (real) ncluster,
1960 &ncluster, TRUE, rlo_bot, rhi_bot);
1964 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1965 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1966 rlo_top, rhi_top, &nlevels);
1969 fprintf(stderr, "\n");
1973 fp = opt2FILE("-om", NFILE, fnm, "w");
1974 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1975 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1976 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1977 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1978 rlo_top, rhi_top, &nlevels);
1983 /* now show what we've done */
1984 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1985 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1986 if (method == m_diagonalize)
1988 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1990 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1993 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1994 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1995 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1997 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);