3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trxio.h"
51 #include "gmx_random.h"
55 #include "gromacs/fileio/futil.h"
56 #include "gromacs/fileio/matio.h"
59 #include "gromacs/fileio/trnio.h"
63 #include "gromacs/linearalgebra/eigensolver.h"
65 /* print to two file pointers at once (i.e. stderr and log) */
67 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
69 fprintf(fp1, "%s", buf);
70 fprintf(fp2, "%s", buf);
73 /* just print a prepared buffer to fp1 and fp2 */
75 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
77 lo_ffprintf(fp1, fp2, buf);
80 /* prepare buffer with one argument, then print to fp1 and fp2 */
82 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
84 sprintf(buf, fmt, arg);
85 lo_ffprintf(fp1, fp2, buf);
88 /* prepare buffer with one argument, then print to fp1 and fp2 */
90 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
92 sprintf(buf, fmt, arg);
93 lo_ffprintf(fp1, fp2, buf);
96 /* prepare buffer with one argument, then print to fp1 and fp2 */
98 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
100 sprintf(buf, fmt, arg);
101 lo_ffprintf(fp1, fp2, buf);
104 /* prepare buffer with two arguments, then print to fp1 and fp2 */
106 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
108 sprintf(buf, fmt, arg1, arg2);
109 lo_ffprintf(fp1, fp2, buf);
112 /* prepare buffer with two arguments, then print to fp1 and fp2 */
114 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
116 sprintf(buf, fmt, arg1, arg2);
117 lo_ffprintf(fp1, fp2, buf);
120 /* prepare buffer with two arguments, then print to fp1 and fp2 */
122 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
124 sprintf(buf, fmt, arg1, arg2);
125 lo_ffprintf(fp1, fp2, buf);
138 void pr_energy(FILE *fp, real e)
140 fprintf(fp, "Energy: %8.4f\n", e);
143 void cp_index(int nn, int from[], int to[])
147 for (i = 0; (i < nn); i++)
153 void mc_optimize(FILE *log, t_mat *m, real *time,
154 int maxiter, int nrandom,
156 const char *conv, output_env_t oenv)
159 real ecur, enext, emin, prob, enorm;
160 int i, j, iswap, jswap, nn, nuphill = 0;
166 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
169 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
170 printf("by reordering the frames to minimize the path between the two structures\n");
171 printf("that have the largest pairwise RMSD.\n");
174 enorm = m->mat[0][0];
175 for (i = 0; (i < m->n1); i++)
177 for (j = 0; (j < m->nn); j++)
179 if (m->mat[i][j] > enorm)
181 enorm = m->mat[i][j];
187 if ((iswap == -1) || (jswap == -1))
189 fprintf(stderr, "Matrix contains identical values in all fields\n");
192 swap_rows(m, 0, iswap);
193 swap_rows(m, m->n1-1, jswap);
194 emin = ecur = mat_energy(m);
195 printf("Largest distance %g between %d and %d. Energy: %g.\n",
196 enorm, iswap, jswap, emin);
198 rng = gmx_rng_init(seed);
201 /* Initiate and store global minimum */
202 minimum = init_mat(nn, m->b1D);
204 copy_t_mat(minimum, m);
208 fp = xvgropen(conv, "Convergence of the MC optimization",
209 "Energy", "Step", oenv);
211 for (i = 0; (i < maxiter); i++)
213 /* Generate new swapping candidates */
216 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
217 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
219 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
221 /* Apply swap and compute energy */
222 swap_rows(m, iswap, jswap);
223 enext = mat_energy(m);
225 /* Compute probability */
227 if ((enext < ecur) || (i < nrandom))
232 /* Store global minimum */
233 copy_t_mat(minimum, m);
239 /* Try Monte Carlo step */
240 prob = exp(-(enext-ecur)/(enorm*kT));
243 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
250 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
251 i, iswap, jswap, enext, prob);
254 fprintf(fp, "%6d %10g\n", i, enext);
260 swap_rows(m, jswap, iswap);
263 fprintf(log, "%d uphill steps were taken during optimization\n",
266 /* Now swap the matrix to get it into global minimum mode */
267 copy_t_mat(m, minimum);
269 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
270 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
271 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
272 for (i = 0; (i < m->nn); i++)
274 fprintf(log, "%10g %5d %10g\n",
277 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
286 static void calc_dist(int nind, rvec x[], real **d)
292 for (i = 0; (i < nind-1); i++)
295 for (j = i+1; (j < nind); j++)
297 /* Should use pbc_dx when analysing multiple molecueles,
298 * but the box is not stored for every frame.
300 rvec_sub(xi, x[j], dx);
306 static real rms_dist(int isize, real **d, real **d_r)
312 for (i = 0; (i < isize-1); i++)
314 for (j = i+1; (j < isize); j++)
316 r = d[i][j]-d_r[i][j];
320 r2 /= (isize*(isize-1))/2;
325 static int rms_dist_comp(const void *a, const void *b)
332 if (da->dist - db->dist < 0)
336 else if (da->dist - db->dist > 0)
343 static int clust_id_comp(const void *a, const void *b)
350 return da->clust - db->clust;
353 static int nrnb_comp(const void *a, const void *b)
360 /* return the b-a, we want highest first */
361 return db->nr - da->nr;
364 void gather(t_mat *m, real cutoff, t_clusters *clust)
368 int i, j, k, nn, cid, n1, diff;
371 /* First we sort the entries in the RMSD matrix */
375 for (i = k = 0; (i < n1); i++)
377 for (j = i+1; (j < n1); j++, k++)
381 d[k].dist = m->mat[i][j];
386 gmx_incons("gather algortihm");
388 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
390 /* Now we make a cluster index for all of the conformations */
393 /* Now we check the closest structures, and equalize their cluster numbers */
394 fprintf(stderr, "Linking structures ");
397 fprintf(stderr, "*");
399 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
401 diff = c[d[k].j].clust - c[d[k].i].clust;
407 c[d[k].j].clust = c[d[k].i].clust;
411 c[d[k].i].clust = c[d[k].j].clust;
417 fprintf(stderr, "\nSorting and renumbering clusters\n");
418 /* Sort on cluster number */
419 qsort(c, n1, sizeof(c[0]), clust_id_comp);
421 /* Renumber clusters */
423 for (k = 1; k < n1; k++)
425 if (c[k].clust != c[k-1].clust)
438 for (k = 0; (k < n1); k++)
440 fprintf(debug, "Cluster index for conformation %d: %d\n",
441 c[k].conf, c[k].clust);
445 for (k = 0; k < n1; k++)
447 clust->cl[c[k].conf] = c[k].clust;
454 gmx_bool jp_same(int **nnb, int i, int j, int P)
460 for (k = 0; nnb[i][k] >= 0; k++)
462 bIn = bIn || (nnb[i][k] == j);
470 for (k = 0; nnb[j][k] >= 0; k++)
472 bIn = bIn || (nnb[j][k] == i);
480 for (ii = 0; nnb[i][ii] >= 0; ii++)
482 for (jj = 0; nnb[j][jj] >= 0; jj++)
484 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
494 static void jarvis_patrick(int n1, real **mat, int M, int P,
495 real rmsdcut, t_clusters *clust)
500 int i, j, k, cid, diff, max;
509 /* First we sort the entries in the RMSD matrix row by row.
510 * This gives us the nearest neighbor list.
514 for (i = 0; (i < n1); i++)
516 for (j = 0; (j < n1); j++)
519 row[j].dist = mat[i][j];
521 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
524 /* Put the M nearest neighbors in the list */
526 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
530 nnb[i][k] = row[j].j;
538 /* Put all neighbors nearer than rmsdcut in the list */
541 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
550 nnb[i][k] = row[j].j;
556 srenew(nnb[i], max+1);
564 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
565 for (i = 0; (i < n1); i++)
567 fprintf(debug, "i:%5d nbs:", i);
568 for (j = 0; nnb[i][j] >= 0; j++)
570 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
572 fprintf(debug, "\n");
577 fprintf(stderr, "Linking structures ");
578 /* Use mcpy for temporary storage of booleans */
579 mcpy = mk_matrix(n1, n1, FALSE);
580 for (i = 0; i < n1; i++)
582 for (j = i+1; j < n1; j++)
584 mcpy[i][j] = jp_same(nnb, i, j, P);
589 fprintf(stderr, "*");
591 for (i = 0; i < n1; i++)
593 for (j = i+1; j < n1; j++)
597 diff = c[j].clust - c[i].clust;
603 c[j].clust = c[i].clust;
607 c[i].clust = c[j].clust;
616 fprintf(stderr, "\nSorting and renumbering clusters\n");
617 /* Sort on cluster number */
618 qsort(c, n1, sizeof(c[0]), clust_id_comp);
620 /* Renumber clusters */
622 for (k = 1; k < n1; k++)
624 if (c[k].clust != c[k-1].clust)
636 for (k = 0; k < n1; k++)
638 clust->cl[c[k].conf] = c[k].clust;
642 for (k = 0; (k < n1); k++)
644 fprintf(debug, "Cluster index for conformation %d: %d\n",
645 c[k].conf, c[k].clust);
649 /* Again, I don't see the point in this... (AF) */
650 /* for(i=0; (i<n1); i++) { */
651 /* for(j=0; (j<n1); j++) */
652 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
654 /* for(i=0; (i<n1); i++) { */
655 /* for(j=0; (j<n1); j++) */
656 /* mat[i][j] = mcpy[i][j]; */
658 done_matrix(n1, &mcpy);
661 for (i = 0; (i < n1); i++)
668 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
672 /* dump neighbor list */
673 fprintf(fp, "%s", title);
674 for (i = 0; (i < n1); i++)
676 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
677 for (j = 0; j < nnb[i].nr; j++)
679 fprintf(fp, "%5d", nnb[i].nb[j]);
685 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
689 int i, j, k, j1, max;
691 /* Put all neighbors nearer than rmsdcut in the list */
692 fprintf(stderr, "Making list of neighbors within cutoff ");
695 for (i = 0; (i < n1); i++)
699 /* put all neighbors within cut-off in list */
700 for (j = 0; j < n1; j++)
702 if (mat[i][j] < rmsdcut)
707 srenew(nnb[i].nb, max);
713 /* store nr of neighbors, we'll need that */
715 if (i%(1+n1/100) == 0)
717 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
720 fprintf(stderr, "%3d%%\n", 100);
723 /* sort neighbor list on number of neighbors, largest first */
724 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
728 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
731 /* turn first structure with all its neighbors (largest) into cluster
732 remove them from pool of structures and repeat for all remaining */
733 fprintf(stderr, "Finding clusters %4d", 0);
734 /* cluster id's start at 1: */
738 /* set cluster id (k) for first item in neighborlist */
739 for (j = 0; j < nnb[0].nr; j++)
741 clust->cl[nnb[0].nb[j]] = k;
747 /* adjust number of neighbors for others, taking removals into account: */
748 for (i = 1; i < n1 && nnb[i].nr; i++)
751 for (j = 0; j < nnb[i].nr; j++)
753 /* if this neighbor wasn't removed */
754 if (clust->cl[nnb[i].nb[j]] == 0)
756 /* shift the rest (j1<=j) */
757 nnb[i].nb[j1] = nnb[i].nb[j];
762 /* now j1 is the new number of neighbors */
765 /* sort again on nnb[].nr, because we have new # neighbors: */
766 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
767 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
769 fprintf(stderr, "\b\b\b\b%4d", k);
773 fprintf(stderr, "\n");
777 fprintf(debug, "Clusters (%d):\n", k);
778 for (i = 0; i < n1; i++)
780 fprintf(debug, " %3d", clust->cl[i]);
782 fprintf(debug, "\n");
788 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
789 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
794 int i, i0, j, max_nf;
802 natom = read_first_x(oenv, &status, fn, &t, &x, box);
809 gmx_rmpbc(gpbc, natom, box, x);
815 srenew(*time, max_nf);
820 /* Store only the interesting atoms */
821 for (j = 0; (j < isize); j++)
823 copy_rvec(x[index[j]], xx[i0][j]);
830 while (read_next_x(oenv, status, &t, x, box));
831 fprintf(stderr, "Allocated %lu bytes for frames\n",
832 (unsigned long) (max_nf*isize*sizeof(**xx)));
833 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
840 static int plot_clusters(int nf, real **mat, t_clusters *clust,
843 int i, j, ncluster, ci;
844 int *cl_id, *nstruct, *strind;
849 for (i = 0; i < nf; i++)
852 cl_id[i] = clust->cl[i];
856 for (i = 0; i < nf; i++)
858 if (nstruct[i] >= minstruct)
861 for (j = 0; (j < nf); j++)
865 strind[j] = ncluster;
871 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
872 ncluster, minstruct);
874 for (i = 0; (i < nf); i++)
877 for (j = 0; j < i; j++)
879 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
881 /* color different clusters with different colors, as long as
882 we don't run out of colors */
883 mat[i][j] = strind[i];
898 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
902 for (i = 0; i < nf; i++)
904 for (j = 0; j < i; j++)
906 if (clust->cl[i] == clust->cl[j])
918 static char *parse_filename(const char *fn, int maxnr)
927 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
929 /* number of digits needed in numbering */
930 i = (int)(log(maxnr)/log(10)) + 1;
931 /* split fn and ext */
932 ext = strrchr(fn, '.');
935 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
938 /* insert e.g. '%03d' between fn and ext */
939 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
940 snew(fnout, strlen(buf)+1);
946 static void ana_trans(t_clusters *clust, int nf,
947 const char *transfn, const char *ntransfn, FILE *log,
948 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
953 int i, ntranst, maxtrans;
956 snew(ntrans, clust->ncl);
957 snew(trans, clust->ncl);
958 snew(axis, clust->ncl);
959 for (i = 0; i < clust->ncl; i++)
962 snew(trans[i], clust->ncl);
966 for (i = 1; i < nf; i++)
968 if (clust->cl[i] != clust->cl[i-1])
971 ntrans[clust->cl[i-1]-1]++;
972 ntrans[clust->cl[i]-1]++;
973 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
974 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
977 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
978 "max %d between two specific clusters\n", ntranst, maxtrans);
981 fp = ffopen(transfn, "w");
982 i = min(maxtrans+1, 80);
983 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
984 "from cluster", "to cluster",
985 clust->ncl, clust->ncl, axis, axis, trans,
986 0, maxtrans, rlo, rhi, &i);
991 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
993 for (i = 0; i < clust->ncl; i++)
995 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1000 for (i = 0; i < clust->ncl; i++)
1008 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1009 int natom, t_atoms *atoms, rvec *xtps,
1010 real *mass, rvec **xx, real *time,
1011 int ifsize, atom_id *fitidx,
1012 int iosize, atom_id *outidx,
1013 const char *trxfn, const char *sizefn,
1014 const char *transfn, const char *ntransfn,
1015 const char *clustidfn, gmx_bool bAverage,
1016 int write_ncl, int write_nst, real rmsmin,
1017 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1018 const output_env_t oenv)
1021 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1022 t_trxstatus *trxout = NULL;
1023 t_trxstatus *trxsout = NULL;
1024 int i, i1, cl, nstr, *structure, first = 0, midstr;
1025 gmx_bool *bWrite = NULL;
1026 real r, clrmsd, midrmsd;
1032 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1036 /* do we write all structures? */
1039 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1042 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1043 bAverage ? "average" : "middle", trxfn);
1046 /* find out what we want to tell the user:
1047 Writing [all structures|structures with rmsd > %g] for
1048 {all|first %d} clusters {with more than %d structures} to %s */
1051 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1055 sprintf(buf1, "all structures");
1057 buf2[0] = buf3[0] = '\0';
1058 if (write_ncl >= clust->ncl)
1062 sprintf(buf2, "all ");
1067 sprintf(buf2, "the first %d ", write_ncl);
1071 sprintf(buf3, " with more than %d structures", write_nst);
1073 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1074 ffprintf(stderr, log, buf);
1077 /* Prepare a reference structure for the orientation of the clusters */
1080 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1082 trxout = open_trx(trxfn, "w");
1083 /* Calculate the average structure in each cluster, *
1084 * all structures are fitted to the first struture of the cluster */
1088 if (transfn || ntransfn)
1090 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1095 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1096 fprintf(fp, "@ s0 symbol 2\n");
1097 fprintf(fp, "@ s0 symbol size 0.2\n");
1098 fprintf(fp, "@ s0 linestyle 0\n");
1099 for (i = 0; i < nf; i++)
1101 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1107 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1108 fprintf(fp, "@g%d type %s\n", 0, "bar");
1110 snew(structure, nf);
1111 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1112 "cl.", "#st", "rmsd", "middle", "rmsd");
1113 for (cl = 1; cl <= clust->ncl; cl++)
1115 /* prepare structures (fit, middle, average) */
1118 for (i = 0; i < natom; i++)
1124 for (i1 = 0; i1 < nf; i1++)
1126 if (clust->cl[i1] == cl)
1128 structure[nstr] = i1;
1130 if (trxfn && (bAverage || write_ncl) )
1134 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1142 do_fit(natom, mass, xx[first], xx[i1]);
1146 for (i = 0; i < natom; i++)
1148 rvec_inc(xav[i], xx[i1][i]);
1156 fprintf(fp, "%8d %8d\n", cl, nstr);
1161 for (i1 = 0; i1 < nstr; i1++)
1166 for (i = 0; i < nstr; i++)
1170 r += rmsd[structure[i]][structure[i1]];
1174 r += rmsd[structure[i1]][structure[i]];
1181 midstr = structure[i1];
1188 /* dump cluster info to logfile */
1191 sprintf(buf1, "%6.3f", clrmsd);
1196 sprintf(buf2, "%5.3f", midrmsd);
1204 sprintf(buf1, "%5s", "");
1205 sprintf(buf2, "%5s", "");
1207 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1208 for (i = 0; i < nstr; i++)
1210 if ((i % 7 == 0) && i)
1212 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1219 fprintf(log, "%s %6g", buf, time[i1]);
1223 /* write structures to trajectory file(s) */
1228 for (i = 0; i < nstr; i++)
1233 if (cl < write_ncl+1 && nstr > write_nst)
1235 /* Dump all structures for this cluster */
1236 /* generate numbered filename (there is a %d in trxfn!) */
1237 sprintf(buf, trxsfn, cl);
1238 trxsout = open_trx(buf, "w");
1239 for (i = 0; i < nstr; i++)
1244 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1248 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1254 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1255 xx[structure[i]], NULL, NULL);
1260 /* Dump the average structure for this cluster */
1263 for (i = 0; i < natom; i++)
1265 svmul(1.0/nstr, xav[i], xav[i]);
1270 for (i = 0; i < natom; i++)
1272 copy_rvec(xx[midstr][i], xav[i]);
1276 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1281 do_fit(natom, mass, xtps, xav);
1284 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1304 static void convert_mat(t_matrix *mat, t_mat *rms)
1309 matrix2real(mat, rms->mat);
1310 /* free input xpm matrix data */
1311 for (i = 0; i < mat->nx; i++)
1313 sfree(mat->matrix[i]);
1317 for (i = 0; i < mat->nx; i++)
1319 for (j = i; j < mat->nx; j++)
1321 rms->sumrms += rms->mat[i][j];
1322 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1325 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1332 int gmx_cluster(int argc, char *argv[])
1334 const char *desc[] = {
1335 "[TT]g_cluster[tt] can cluster structures using several different methods.",
1336 "Distances between structures can be determined from a trajectory",
1337 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1338 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1339 "can be used to define the distance between structures.[PAR]",
1341 "single linkage: add a structure to a cluster when its distance to any",
1342 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1344 "Jarvis Patrick: add a structure to a cluster when this structure",
1345 "and a structure in the cluster have each other as neighbors and",
1346 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1347 "of a structure are the M closest structures or all structures within",
1348 "[TT]cutoff[tt].[PAR]",
1350 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1351 "the order of the frames is using the smallest possible increments.",
1352 "With this it is possible to make a smooth animation going from one",
1353 "structure to another with the largest possible (e.g.) RMSD between",
1354 "them, however the intermediate steps should be as small as possible.",
1355 "Applications could be to visualize a potential of mean force",
1356 "ensemble of simulations or a pulling simulation. Obviously the user",
1357 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1358 "The final result can be inspect visually by looking at the matrix",
1359 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1361 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1363 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1364 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1365 "Count number of neighbors using cut-off, take structure with",
1366 "largest number of neighbors with all its neighbors as cluster",
1367 "and eliminate it from the pool of clusters. Repeat for remaining",
1368 "structures in pool.[PAR]",
1370 "When the clustering algorithm assigns each structure to exactly one",
1371 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1372 "file is supplied, the structure with",
1373 "the smallest average distance to the others or the average structure",
1374 "or all structures for each cluster will be written to a trajectory",
1375 "file. When writing all structures, separate numbered files are made",
1376 "for each cluster.[PAR]",
1378 "Two output files are always written:[BR]",
1379 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1380 "and a graphical depiction of the clusters in the lower right half",
1381 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1382 "when two structures are in the same cluster.",
1383 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1385 "[TT]-g[tt] writes information on the options used and a detailed list",
1386 "of all clusters and their members.[PAR]",
1388 "Additionally, a number of optional output files can be written:[BR]",
1389 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1390 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1391 "diagonalization.[BR]",
1392 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1393 "[TT]-tr[tt] writes a matrix of the number transitions between",
1394 "cluster pairs.[BR]",
1395 "[TT]-ntr[tt] writes the total number of transitions to or from",
1396 "each cluster.[BR]",
1397 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1398 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1399 "structure of each cluster or writes numbered files with cluster members",
1400 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1401 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1402 "structure with the smallest average RMSD from all other structures",
1403 "of the cluster.[BR]",
1407 int nf, i, i1, i2, j;
1408 gmx_large_int_t nrms = 0;
1411 rvec *xtps, *usextps, *x1, **xx = NULL;
1412 const char *fn, *trx_out_fn;
1414 t_mat *rms, *orig = NULL;
1419 t_matrix *readmat = NULL;
1422 int isize = 0, ifsize = 0, iosize = 0;
1423 atom_id *index = NULL, *fitidx, *outidx;
1425 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1426 char buf[STRLEN], buf1[80], title[STRLEN];
1427 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1429 int method, ncluster = 0;
1430 static const char *methodname[] = {
1431 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1432 "diagonalization", "gromos", NULL
1435 m_null, m_linkage, m_jarvis_patrick,
1436 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1438 /* Set colors for plotting: white = zero RMS, black = maximum */
1439 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1440 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1441 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1442 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1443 static int nlevels = 40, skip = 1;
1444 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1445 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1446 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1447 static real kT = 1e-3;
1448 static int M = 10, P = 3;
1450 gmx_rmpbc_t gpbc = NULL;
1453 { "-dista", FALSE, etBOOL, {&bRMSdist},
1454 "Use RMSD of distances instead of RMS deviation" },
1455 { "-nlevels", FALSE, etINT, {&nlevels},
1456 "Discretize RMSD matrix in this number of levels" },
1457 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1458 "RMSD cut-off (nm) for two structures to be neighbor" },
1459 { "-fit", FALSE, etBOOL, {&bFit},
1460 "Use least squares fitting before RMSD calculation" },
1461 { "-max", FALSE, etREAL, {&scalemax},
1462 "Maximum level in RMSD matrix" },
1463 { "-skip", FALSE, etINT, {&skip},
1464 "Only analyze every nr-th frame" },
1465 { "-av", FALSE, etBOOL, {&bAverage},
1466 "Write average iso middle structure for each cluster" },
1467 { "-wcl", FALSE, etINT, {&write_ncl},
1468 "Write the structures for this number of clusters to numbered files" },
1469 { "-nst", FALSE, etINT, {&write_nst},
1470 "Only write all structures if more than this number of structures per cluster" },
1471 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1472 "minimum rms difference with rest of cluster for writing structures" },
1473 { "-method", FALSE, etENUM, {methodname},
1474 "Method for cluster determination" },
1475 { "-minstruct", FALSE, etINT, {&minstruct},
1476 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1477 { "-binary", FALSE, etBOOL, {&bBinary},
1478 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1479 "is given by [TT]-cutoff[tt]" },
1480 { "-M", FALSE, etINT, {&M},
1481 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1482 "0 is use cutoff" },
1483 { "-P", FALSE, etINT, {&P},
1484 "Number of identical nearest neighbors required to form a cluster" },
1485 { "-seed", FALSE, etINT, {&seed},
1486 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1487 { "-niter", FALSE, etINT, {&niter},
1488 "Number of iterations for MC" },
1489 { "-nrandom", FALSE, etINT, {&nrandom},
1490 "The first iterations for MC may be done complete random, to shuffle the frames" },
1491 { "-kT", FALSE, etREAL, {&kT},
1492 "Boltzmann weighting factor for Monte Carlo optimization "
1493 "(zero turns off uphill steps)" },
1494 { "-pbc", FALSE, etBOOL,
1495 { &bPBC }, "PBC check" }
1498 { efTRX, "-f", NULL, ffOPTRD },
1499 { efTPS, "-s", NULL, ffOPTRD },
1500 { efNDX, NULL, NULL, ffOPTRD },
1501 { efXPM, "-dm", "rmsd", ffOPTRD },
1502 { efXPM, "-om", "rmsd-raw", ffWRITE },
1503 { efXPM, "-o", "rmsd-clust", ffWRITE },
1504 { efLOG, "-g", "cluster", ffWRITE },
1505 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1506 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1507 { efXVG, "-conv", "mc-conv", ffOPTWR },
1508 { efXVG, "-sz", "clust-size", ffOPTWR},
1509 { efXPM, "-tr", "clust-trans", ffOPTWR},
1510 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1511 { efXVG, "-clid", "clust-id.xvg", ffOPTWR},
1512 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1514 #define NFILE asize(fnm)
1516 if (!parse_common_args(&argc, argv,
1517 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1518 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1525 bReadMat = opt2bSet("-dm", NFILE, fnm);
1526 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1527 if (opt2parg_bSet("-av", asize(pa), pa) ||
1528 opt2parg_bSet("-wcl", asize(pa), pa) ||
1529 opt2parg_bSet("-nst", asize(pa), pa) ||
1530 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1531 opt2bSet("-cl", NFILE, fnm) )
1533 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1539 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1542 "\nWarning: assuming the time unit in %s is %s\n",
1543 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1545 if (trx_out_fn && !bReadTraj)
1547 fprintf(stderr, "\nWarning: "
1548 "cannot write cluster structures without reading trajectory\n"
1549 " ignoring option -cl %s\n", trx_out_fn);
1553 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1559 gmx_fatal(FARGS, "Invalid method");
1562 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1563 method == m_gromos );
1566 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1568 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1569 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1571 /* check input and write parameters to log file */
1572 bUseRmsdCut = FALSE;
1573 if (method == m_jarvis_patrick)
1575 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1576 if ((M < 0) || (M == 1))
1578 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1582 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1589 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1593 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1598 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1601 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1603 else /* method != m_jarvis */
1605 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1607 if (bUseRmsdCut && method != m_jarvis_patrick)
1609 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1611 if (method == m_monte_carlo)
1613 fprintf(log, "Using %d iterations\n", niter);
1618 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1624 /* don't read mass-database as masses (and top) are not used */
1625 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1629 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1632 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1633 bReadMat ? "" : " and RMSD calculation");
1634 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1635 1, &ifsize, &fitidx, &grpname);
1638 fprintf(stderr, "\nSelect group for output:\n");
1639 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1640 1, &iosize, &outidx, &grpname);
1641 /* merge and convert both index groups: */
1642 /* first copy outidx to index. let outidx refer to elements in index */
1643 snew(index, iosize);
1645 for (i = 0; i < iosize; i++)
1647 index[i] = outidx[i];
1650 /* now lookup elements from fitidx in index, add them if necessary
1651 and also let fitidx refer to elements in index */
1652 for (i = 0; i < ifsize; i++)
1655 while (j < isize && index[j] != fitidx[i])
1661 /* slow this way, but doesn't matter much */
1663 srenew(index, isize);
1665 index[j] = fitidx[i];
1669 else /* !trx_out_fn */
1673 for (i = 0; i < ifsize; i++)
1675 index[i] = fitidx[i];
1683 /* Loop over first coordinate file */
1684 fn = opt2fn("-f", NFILE, fnm);
1686 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1687 output_env_conv_times(oenv, nf, time);
1688 if (!bRMSdist || bAnalyze)
1690 /* Center all frames on zero */
1692 for (i = 0; i < ifsize; i++)
1694 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1698 for (i = 0; i < nf; i++)
1700 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1706 gmx_rmpbc_done(gpbc);
1712 fprintf(stderr, "Reading rms distance matrix ");
1713 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1714 fprintf(stderr, "\n");
1715 if (readmat[0].nx != readmat[0].ny)
1717 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1718 readmat[0].nx, readmat[0].ny);
1720 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1722 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1723 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1728 time = readmat[0].axis_x;
1729 time_invfac = output_env_get_time_invfactor(oenv);
1730 for (i = 0; i < nf; i++)
1732 time[i] *= time_invfac;
1735 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1736 convert_mat(&(readmat[0]), rms);
1738 nlevels = readmat[0].nmap;
1740 else /* !bReadMat */
1742 rms = init_mat(nf, method == m_diagonalize);
1743 nrms = ((gmx_large_int_t)nf*((gmx_large_int_t)nf-1))/2;
1746 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1747 /* Initialize work array */
1749 for (i1 = 0; i1 < nf; i1++)
1751 for (i2 = i1+1; i2 < nf; i2++)
1753 for (i = 0; i < isize; i++)
1755 copy_rvec(xx[i1][i], x1[i]);
1759 do_fit(isize, mass, xx[i2], x1);
1761 rmsd = rmsdev(isize, mass, xx[i2], x1);
1762 set_mat_entry(rms, i1, i2, rmsd);
1764 nrms -= (gmx_large_int_t) (nf-i1-1);
1765 fprintf(stderr, "\r# RMSD calculations left: "gmx_large_int_pfmt" ", nrms);
1771 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1773 /* Initiate work arrays */
1776 for (i = 0; (i < isize); i++)
1781 for (i1 = 0; i1 < nf ; i1++)
1783 calc_dist(isize, xx[i1], d1);
1784 for (i2 = i1+1; (i2 < nf); i2++)
1786 calc_dist(isize, xx[i2], d2);
1787 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1790 fprintf(stderr, "\r# RMSD calculations left: "gmx_large_int_pfmt" ", nrms);
1792 /* Clean up work arrays */
1793 for (i = 0; (i < isize); i++)
1801 fprintf(stderr, "\n\n");
1803 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1804 rms->minrms, rms->maxrms);
1805 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1806 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1807 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1808 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1810 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1811 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1813 if (bAnalyze && (rmsmin < rms->minrms) )
1815 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1816 rmsmin, rms->minrms);
1818 if (bAnalyze && (rmsmin > rmsdcut) )
1820 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1824 /* Plot the rmsd distribution */
1825 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1829 for (i1 = 0; (i1 < nf); i1++)
1831 for (i2 = 0; (i2 < nf); i2++)
1833 if (rms->mat[i1][i2] < rmsdcut)
1835 rms->mat[i1][i2] = 0;
1839 rms->mat[i1][i2] = 1;
1849 /* Now sort the matrix and write it out again */
1850 gather(rms, rmsdcut, &clust);
1853 /* Do a diagonalization */
1854 snew(eigenvalues, nf);
1855 snew(eigenvectors, nf*nf);
1856 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1857 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1858 sfree(eigenvectors);
1860 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1861 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1862 for (i = 0; (i < nf); i++)
1864 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1869 orig = init_mat(rms->nn, FALSE);
1871 copy_t_mat(orig, rms);
1872 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1873 opt2fn_null("-conv", NFILE, fnm), oenv);
1875 case m_jarvis_patrick:
1876 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1879 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1882 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1885 if (method == m_monte_carlo || method == m_diagonalize)
1887 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1895 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1899 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1901 init_t_atoms(&useatoms, isize, FALSE);
1902 snew(usextps, isize);
1903 useatoms.resinfo = top.atoms.resinfo;
1904 for (i = 0; i < isize; i++)
1906 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1907 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1908 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1909 copy_rvec(xtps[index[i]], usextps[i]);
1911 useatoms.nr = isize;
1912 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1913 ifsize, fitidx, iosize, outidx,
1914 bReadTraj ? trx_out_fn : NULL,
1915 opt2fn_null("-sz", NFILE, fnm),
1916 opt2fn_null("-tr", NFILE, fnm),
1917 opt2fn_null("-ntr", NFILE, fnm),
1918 opt2fn_null("-clid", NFILE, fnm),
1919 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1920 rlo_bot, rhi_bot, oenv);
1924 if (bBinary && !bAnalyze)
1926 /* Make the clustering visible */
1927 for (i2 = 0; (i2 < nf); i2++)
1929 for (i1 = i2+1; (i1 < nf); i1++)
1931 if (rms->mat[i1][i2])
1933 rms->mat[i1][i2] = rms->maxrms;
1939 fp = opt2FILE("-o", NFILE, fnm, "w");
1940 fprintf(stderr, "Writing rms distance/clustering matrix ");
1943 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1944 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1945 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1949 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1950 sprintf(title, "RMS%sDeviation / Cluster Index",
1951 bRMSdist ? " Distance " : " ");
1954 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1955 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1956 rlo_top, rhi_top, 0.0, (real) ncluster,
1957 &ncluster, TRUE, rlo_bot, rhi_bot);
1961 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1962 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1963 rlo_top, rhi_top, &nlevels);
1966 fprintf(stderr, "\n");
1970 fp = opt2FILE("-om", NFILE, fnm, "w");
1971 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1972 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1973 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1974 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1975 rlo_top, rhi_top, &nlevels);
1980 /* now show what we've done */
1981 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1982 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1983 if (method == m_diagonalize)
1985 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1987 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1990 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1991 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1992 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1994 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);