2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/tpxio.h"
46 #include "gromacs/fileio/trnio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/cmat.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/legacyheaders/macros.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/legacyheaders/viewit.h"
54 #include "gromacs/linearalgebra/eigensolver.h"
55 #include "gromacs/math/do_fit.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/pbcutil/pbc.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/random/random.h"
60 #include "gromacs/topology/index.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/smalloc.h"
66 /* print to two file pointers at once (i.e. stderr and log) */
68 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
70 fprintf(fp1, "%s", buf);
71 fprintf(fp2, "%s", buf);
74 /* just print a prepared buffer to fp1 and fp2 */
76 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
78 lo_ffprintf(fp1, fp2, buf);
81 /* prepare buffer with one argument, then print to fp1 and fp2 */
83 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
85 sprintf(buf, fmt, arg);
86 lo_ffprintf(fp1, fp2, buf);
89 /* prepare buffer with one argument, then print to fp1 and fp2 */
91 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
93 sprintf(buf, fmt, arg);
94 lo_ffprintf(fp1, fp2, buf);
97 /* prepare buffer with one argument, then print to fp1 and fp2 */
99 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
101 sprintf(buf, fmt, arg);
102 lo_ffprintf(fp1, fp2, buf);
105 /* prepare buffer with two arguments, then print to fp1 and fp2 */
107 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
109 sprintf(buf, fmt, arg1, arg2);
110 lo_ffprintf(fp1, fp2, buf);
113 /* prepare buffer with two arguments, then print to fp1 and fp2 */
115 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
117 sprintf(buf, fmt, arg1, arg2);
118 lo_ffprintf(fp1, fp2, buf);
121 /* prepare buffer with two arguments, then print to fp1 and fp2 */
123 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
125 sprintf(buf, fmt, arg1, arg2);
126 lo_ffprintf(fp1, fp2, buf);
139 void cp_index(int nn, int from[], int to[])
143 for (i = 0; (i < nn); i++)
149 void mc_optimize(FILE *log, t_mat *m, real *time,
150 int maxiter, int nrandom,
152 const char *conv, output_env_t oenv)
155 real ecur, enext, emin, prob, enorm;
156 int i, j, iswap, jswap, nn, nuphill = 0;
162 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
165 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
166 printf("by reordering the frames to minimize the path between the two structures\n");
167 printf("that have the largest pairwise RMSD.\n");
170 enorm = m->mat[0][0];
171 for (i = 0; (i < m->n1); i++)
173 for (j = 0; (j < m->nn); j++)
175 if (m->mat[i][j] > enorm)
177 enorm = m->mat[i][j];
183 if ((iswap == -1) || (jswap == -1))
185 fprintf(stderr, "Matrix contains identical values in all fields\n");
188 swap_rows(m, 0, iswap);
189 swap_rows(m, m->n1-1, jswap);
190 emin = ecur = mat_energy(m);
191 printf("Largest distance %g between %d and %d. Energy: %g.\n",
192 enorm, iswap, jswap, emin);
194 rng = gmx_rng_init(seed);
197 /* Initiate and store global minimum */
198 minimum = init_mat(nn, m->b1D);
200 copy_t_mat(minimum, m);
204 fp = xvgropen(conv, "Convergence of the MC optimization",
205 "Energy", "Step", oenv);
207 for (i = 0; (i < maxiter); i++)
209 /* Generate new swapping candidates */
212 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
213 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
215 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
217 /* Apply swap and compute energy */
218 swap_rows(m, iswap, jswap);
219 enext = mat_energy(m);
221 /* Compute probability */
223 if ((enext < ecur) || (i < nrandom))
228 /* Store global minimum */
229 copy_t_mat(minimum, m);
235 /* Try Monte Carlo step */
236 prob = exp(-(enext-ecur)/(enorm*kT));
239 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
246 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
247 i, iswap, jswap, enext, prob);
250 fprintf(fp, "%6d %10g\n", i, enext);
256 swap_rows(m, jswap, iswap);
259 fprintf(log, "%d uphill steps were taken during optimization\n",
262 /* Now swap the matrix to get it into global minimum mode */
263 copy_t_mat(m, minimum);
265 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
266 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
267 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
268 for (i = 0; (i < m->nn); i++)
270 fprintf(log, "%10g %5d %10g\n",
273 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
282 static void calc_dist(int nind, rvec x[], real **d)
288 for (i = 0; (i < nind-1); i++)
291 for (j = i+1; (j < nind); j++)
293 /* Should use pbc_dx when analysing multiple molecueles,
294 * but the box is not stored for every frame.
296 rvec_sub(xi, x[j], dx);
302 static real rms_dist(int isize, real **d, real **d_r)
308 for (i = 0; (i < isize-1); i++)
310 for (j = i+1; (j < isize); j++)
312 r = d[i][j]-d_r[i][j];
316 r2 /= (isize*(isize-1))/2;
321 static int rms_dist_comp(const void *a, const void *b)
328 if (da->dist - db->dist < 0)
332 else if (da->dist - db->dist > 0)
339 static int clust_id_comp(const void *a, const void *b)
346 return da->clust - db->clust;
349 static int nrnb_comp(const void *a, const void *b)
356 /* return the b-a, we want highest first */
357 return db->nr - da->nr;
360 void gather(t_mat *m, real cutoff, t_clusters *clust)
364 int i, j, k, nn, cid, n1, diff;
367 /* First we sort the entries in the RMSD matrix */
371 for (i = k = 0; (i < n1); i++)
373 for (j = i+1; (j < n1); j++, k++)
377 d[k].dist = m->mat[i][j];
382 gmx_incons("gather algortihm");
384 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
386 /* Now we make a cluster index for all of the conformations */
389 /* Now we check the closest structures, and equalize their cluster numbers */
390 fprintf(stderr, "Linking structures ");
393 fprintf(stderr, "*");
395 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
397 diff = c[d[k].j].clust - c[d[k].i].clust;
403 c[d[k].j].clust = c[d[k].i].clust;
407 c[d[k].i].clust = c[d[k].j].clust;
413 fprintf(stderr, "\nSorting and renumbering clusters\n");
414 /* Sort on cluster number */
415 qsort(c, n1, sizeof(c[0]), clust_id_comp);
417 /* Renumber clusters */
419 for (k = 1; k < n1; k++)
421 if (c[k].clust != c[k-1].clust)
434 for (k = 0; (k < n1); k++)
436 fprintf(debug, "Cluster index for conformation %d: %d\n",
437 c[k].conf, c[k].clust);
441 for (k = 0; k < n1; k++)
443 clust->cl[c[k].conf] = c[k].clust;
450 gmx_bool jp_same(int **nnb, int i, int j, int P)
456 for (k = 0; nnb[i][k] >= 0; k++)
458 bIn = bIn || (nnb[i][k] == j);
466 for (k = 0; nnb[j][k] >= 0; k++)
468 bIn = bIn || (nnb[j][k] == i);
476 for (ii = 0; nnb[i][ii] >= 0; ii++)
478 for (jj = 0; nnb[j][jj] >= 0; jj++)
480 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
490 static void jarvis_patrick(int n1, real **mat, int M, int P,
491 real rmsdcut, t_clusters *clust)
496 int i, j, k, cid, diff, max;
505 /* First we sort the entries in the RMSD matrix row by row.
506 * This gives us the nearest neighbor list.
510 for (i = 0; (i < n1); i++)
512 for (j = 0; (j < n1); j++)
515 row[j].dist = mat[i][j];
517 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
520 /* Put the M nearest neighbors in the list */
522 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
526 nnb[i][k] = row[j].j;
534 /* Put all neighbors nearer than rmsdcut in the list */
537 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
546 nnb[i][k] = row[j].j;
552 srenew(nnb[i], max+1);
560 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
561 for (i = 0; (i < n1); i++)
563 fprintf(debug, "i:%5d nbs:", i);
564 for (j = 0; nnb[i][j] >= 0; j++)
566 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
568 fprintf(debug, "\n");
573 fprintf(stderr, "Linking structures ");
574 /* Use mcpy for temporary storage of booleans */
575 mcpy = mk_matrix(n1, n1, FALSE);
576 for (i = 0; i < n1; i++)
578 for (j = i+1; j < n1; j++)
580 mcpy[i][j] = jp_same(nnb, i, j, P);
585 fprintf(stderr, "*");
587 for (i = 0; i < n1; i++)
589 for (j = i+1; j < n1; j++)
593 diff = c[j].clust - c[i].clust;
599 c[j].clust = c[i].clust;
603 c[i].clust = c[j].clust;
612 fprintf(stderr, "\nSorting and renumbering clusters\n");
613 /* Sort on cluster number */
614 qsort(c, n1, sizeof(c[0]), clust_id_comp);
616 /* Renumber clusters */
618 for (k = 1; k < n1; k++)
620 if (c[k].clust != c[k-1].clust)
632 for (k = 0; k < n1; k++)
634 clust->cl[c[k].conf] = c[k].clust;
638 for (k = 0; (k < n1); k++)
640 fprintf(debug, "Cluster index for conformation %d: %d\n",
641 c[k].conf, c[k].clust);
645 /* Again, I don't see the point in this... (AF) */
646 /* for(i=0; (i<n1); i++) { */
647 /* for(j=0; (j<n1); j++) */
648 /* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
650 /* for(i=0; (i<n1); i++) { */
651 /* for(j=0; (j<n1); j++) */
652 /* mat[i][j] = mcpy[i][j]; */
654 done_matrix(n1, &mcpy);
657 for (i = 0; (i < n1); i++)
664 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
668 /* dump neighbor list */
669 fprintf(fp, "%s", title);
670 for (i = 0; (i < n1); i++)
672 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
673 for (j = 0; j < nnb[i].nr; j++)
675 fprintf(fp, "%5d", nnb[i].nb[j]);
681 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
685 int i, j, k, j1, max;
687 /* Put all neighbors nearer than rmsdcut in the list */
688 fprintf(stderr, "Making list of neighbors within cutoff ");
691 for (i = 0; (i < n1); i++)
695 /* put all neighbors within cut-off in list */
696 for (j = 0; j < n1; j++)
698 if (mat[i][j] < rmsdcut)
703 srenew(nnb[i].nb, max);
709 /* store nr of neighbors, we'll need that */
711 if (i%(1+n1/100) == 0)
713 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
716 fprintf(stderr, "%3d%%\n", 100);
719 /* sort neighbor list on number of neighbors, largest first */
720 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
724 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
727 /* turn first structure with all its neighbors (largest) into cluster
728 remove them from pool of structures and repeat for all remaining */
729 fprintf(stderr, "Finding clusters %4d", 0);
730 /* cluster id's start at 1: */
734 /* set cluster id (k) for first item in neighborlist */
735 for (j = 0; j < nnb[0].nr; j++)
737 clust->cl[nnb[0].nb[j]] = k;
743 /* adjust number of neighbors for others, taking removals into account: */
744 for (i = 1; i < n1 && nnb[i].nr; i++)
747 for (j = 0; j < nnb[i].nr; j++)
749 /* if this neighbor wasn't removed */
750 if (clust->cl[nnb[i].nb[j]] == 0)
752 /* shift the rest (j1<=j) */
753 nnb[i].nb[j1] = nnb[i].nb[j];
758 /* now j1 is the new number of neighbors */
761 /* sort again on nnb[].nr, because we have new # neighbors: */
762 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
763 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
765 fprintf(stderr, "\b\b\b\b%4d", k);
769 fprintf(stderr, "\n");
773 fprintf(debug, "Clusters (%d):\n", k);
774 for (i = 0; i < n1; i++)
776 fprintf(debug, " %3d", clust->cl[i]);
778 fprintf(debug, "\n");
784 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
785 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
790 int i, i0, j, max_nf;
798 natom = read_first_x(oenv, &status, fn, &t, &x, box);
805 gmx_rmpbc(gpbc, natom, box, x);
811 srenew(*time, max_nf);
816 /* Store only the interesting atoms */
817 for (j = 0; (j < isize); j++)
819 copy_rvec(x[index[j]], xx[i0][j]);
826 while (read_next_x(oenv, status, &t, x, box));
827 fprintf(stderr, "Allocated %lu bytes for frames\n",
828 (unsigned long) (max_nf*isize*sizeof(**xx)));
829 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
836 static int plot_clusters(int nf, real **mat, t_clusters *clust,
839 int i, j, ncluster, ci;
840 int *cl_id, *nstruct, *strind;
845 for (i = 0; i < nf; i++)
848 cl_id[i] = clust->cl[i];
852 for (i = 0; i < nf; i++)
854 if (nstruct[i] >= minstruct)
857 for (j = 0; (j < nf); j++)
861 strind[j] = ncluster;
867 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
868 ncluster, minstruct);
870 for (i = 0; (i < nf); i++)
873 for (j = 0; j < i; j++)
875 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
877 /* color different clusters with different colors, as long as
878 we don't run out of colors */
879 mat[i][j] = strind[i];
894 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
898 for (i = 0; i < nf; i++)
900 for (j = 0; j < i; j++)
902 if (clust->cl[i] == clust->cl[j])
914 static char *parse_filename(const char *fn, int maxnr)
923 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
925 /* number of digits needed in numbering */
926 i = (int)(log(maxnr)/log(10)) + 1;
927 /* split fn and ext */
928 ext = strrchr(fn, '.');
931 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
934 /* insert e.g. '%03d' between fn and ext */
935 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
936 snew(fnout, strlen(buf)+1);
942 static void ana_trans(t_clusters *clust, int nf,
943 const char *transfn, const char *ntransfn, FILE *log,
944 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
949 int i, ntranst, maxtrans;
952 snew(ntrans, clust->ncl);
953 snew(trans, clust->ncl);
954 snew(axis, clust->ncl);
955 for (i = 0; i < clust->ncl; i++)
958 snew(trans[i], clust->ncl);
962 for (i = 1; i < nf; i++)
964 if (clust->cl[i] != clust->cl[i-1])
967 ntrans[clust->cl[i-1]-1]++;
968 ntrans[clust->cl[i]-1]++;
969 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
970 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
973 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
974 "max %d between two specific clusters\n", ntranst, maxtrans);
977 fp = gmx_ffopen(transfn, "w");
978 i = min(maxtrans+1, 80);
979 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
980 "from cluster", "to cluster",
981 clust->ncl, clust->ncl, axis, axis, trans,
982 0, maxtrans, rlo, rhi, &i);
987 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
989 for (i = 0; i < clust->ncl; i++)
991 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
996 for (i = 0; i < clust->ncl; i++)
1004 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1005 int natom, t_atoms *atoms, rvec *xtps,
1006 real *mass, rvec **xx, real *time,
1007 int ifsize, atom_id *fitidx,
1008 int iosize, atom_id *outidx,
1009 const char *trxfn, const char *sizefn,
1010 const char *transfn, const char *ntransfn,
1011 const char *clustidfn, gmx_bool bAverage,
1012 int write_ncl, int write_nst, real rmsmin,
1013 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1014 const output_env_t oenv)
1016 FILE *size_fp = NULL;
1017 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1018 t_trxstatus *trxout = NULL;
1019 t_trxstatus *trxsout = NULL;
1020 int i, i1, cl, nstr, *structure, first = 0, midstr;
1021 gmx_bool *bWrite = NULL;
1022 real r, clrmsd, midrmsd;
1028 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1032 /* do we write all structures? */
1035 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl));
1038 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1039 bAverage ? "average" : "middle", trxfn);
1042 /* find out what we want to tell the user:
1043 Writing [all structures|structures with rmsd > %g] for
1044 {all|first %d} clusters {with more than %d structures} to %s */
1047 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1051 sprintf(buf1, "all structures");
1053 buf2[0] = buf3[0] = '\0';
1054 if (write_ncl >= clust->ncl)
1058 sprintf(buf2, "all ");
1063 sprintf(buf2, "the first %d ", write_ncl);
1067 sprintf(buf3, " with more than %d structures", write_nst);
1069 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1070 ffprintf(stderr, log, buf);
1073 /* Prepare a reference structure for the orientation of the clusters */
1076 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1078 trxout = open_trx(trxfn, "w");
1079 /* Calculate the average structure in each cluster, *
1080 * all structures are fitted to the first struture of the cluster */
1084 if (transfn || ntransfn)
1086 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1091 FILE *fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1092 if (output_env_get_print_xvgr_codes(oenv))
1094 fprintf(fp, "@ s0 symbol 2\n");
1095 fprintf(fp, "@ s0 symbol size 0.2\n");
1096 fprintf(fp, "@ s0 linestyle 0\n");
1098 for (i = 0; i < nf; i++)
1100 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1106 size_fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1107 if (output_env_get_print_xvgr_codes(oenv))
1109 fprintf(size_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(size_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);
1311 static void convert_mat(t_matrix *mat, t_mat *rms)
1316 matrix2real(mat, rms->mat);
1317 /* free input xpm matrix data */
1318 for (i = 0; i < mat->nx; i++)
1320 sfree(mat->matrix[i]);
1324 for (i = 0; i < mat->nx; i++)
1326 for (j = i; j < mat->nx; j++)
1328 rms->sumrms += rms->mat[i][j];
1329 rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
1332 rms->minrms = min(rms->minrms, rms->mat[i][j]);
1339 int gmx_cluster(int argc, char *argv[])
1341 const char *desc[] = {
1342 "[THISMODULE] can cluster structures using several different methods.",
1343 "Distances between structures can be determined from a trajectory",
1344 "or read from an [REF].xpm[ref] matrix file with the [TT]-dm[tt] option.",
1345 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1346 "can be used to define the distance between structures.[PAR]",
1348 "single linkage: add a structure to a cluster when its distance to any",
1349 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1351 "Jarvis Patrick: add a structure to a cluster when this structure",
1352 "and a structure in the cluster have each other as neighbors and",
1353 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1354 "of a structure are the M closest structures or all structures within",
1355 "[TT]cutoff[tt].[PAR]",
1357 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1358 "the order of the frames is using the smallest possible increments.",
1359 "With this it is possible to make a smooth animation going from one",
1360 "structure to another with the largest possible (e.g.) RMSD between",
1361 "them, however the intermediate steps should be as small as possible.",
1362 "Applications could be to visualize a potential of mean force",
1363 "ensemble of simulations or a pulling simulation. Obviously the user",
1364 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1365 "The final result can be inspect visually by looking at the matrix",
1366 "[REF].xpm[ref] file, which should vary smoothly from bottom to top.[PAR]",
1368 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1370 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1371 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1372 "Count number of neighbors using cut-off, take structure with",
1373 "largest number of neighbors with all its neighbors as cluster",
1374 "and eliminate it from the pool of clusters. Repeat for remaining",
1375 "structures in pool.[PAR]",
1377 "When the clustering algorithm assigns each structure to exactly one",
1378 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1379 "file is supplied, the structure with",
1380 "the smallest average distance to the others or the average structure",
1381 "or all structures for each cluster will be written to a trajectory",
1382 "file. When writing all structures, separate numbered files are made",
1383 "for each cluster.[PAR]",
1385 "Two output files are always written:[BR]",
1386 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1387 "and a graphical depiction of the clusters in the lower right half",
1388 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1389 "when two structures are in the same cluster.",
1390 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1392 "[TT]-g[tt] writes information on the options used and a detailed list",
1393 "of all clusters and their members.[PAR]",
1395 "Additionally, a number of optional output files can be written:[BR]",
1396 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1397 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1398 "diagonalization.[BR]",
1399 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1400 "[TT]-tr[tt] writes a matrix of the number transitions between",
1401 "cluster pairs.[BR]",
1402 "[TT]-ntr[tt] writes the total number of transitions to or from",
1403 "each cluster.[BR]",
1404 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1405 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1406 "structure of each cluster or writes numbered files with cluster members",
1407 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1408 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1409 "structure with the smallest average RMSD from all other structures",
1410 "of the cluster.[BR]",
1414 int nf, i, i1, i2, j;
1415 gmx_int64_t nrms = 0;
1418 rvec *xtps, *usextps, *x1, **xx = NULL;
1419 const char *fn, *trx_out_fn;
1421 t_mat *rms, *orig = NULL;
1426 t_matrix *readmat = NULL;
1429 int isize = 0, ifsize = 0, iosize = 0;
1430 atom_id *index = NULL, *fitidx, *outidx;
1432 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1433 char buf[STRLEN], buf1[80], title[STRLEN];
1434 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1436 int method, ncluster = 0;
1437 static const char *methodname[] = {
1438 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1439 "diagonalization", "gromos", NULL
1442 m_null, m_linkage, m_jarvis_patrick,
1443 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1445 /* Set colors for plotting: white = zero RMS, black = maximum */
1446 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1447 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1448 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1449 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1450 static int nlevels = 40, skip = 1;
1451 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1452 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1453 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1454 static real kT = 1e-3;
1455 static int M = 10, P = 3;
1457 gmx_rmpbc_t gpbc = NULL;
1460 { "-dista", FALSE, etBOOL, {&bRMSdist},
1461 "Use RMSD of distances instead of RMS deviation" },
1462 { "-nlevels", FALSE, etINT, {&nlevels},
1463 "Discretize RMSD matrix in this number of levels" },
1464 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1465 "RMSD cut-off (nm) for two structures to be neighbor" },
1466 { "-fit", FALSE, etBOOL, {&bFit},
1467 "Use least squares fitting before RMSD calculation" },
1468 { "-max", FALSE, etREAL, {&scalemax},
1469 "Maximum level in RMSD matrix" },
1470 { "-skip", FALSE, etINT, {&skip},
1471 "Only analyze every nr-th frame" },
1472 { "-av", FALSE, etBOOL, {&bAverage},
1473 "Write average iso middle structure for each cluster" },
1474 { "-wcl", FALSE, etINT, {&write_ncl},
1475 "Write the structures for this number of clusters to numbered files" },
1476 { "-nst", FALSE, etINT, {&write_nst},
1477 "Only write all structures if more than this number of structures per cluster" },
1478 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1479 "minimum rms difference with rest of cluster for writing structures" },
1480 { "-method", FALSE, etENUM, {methodname},
1481 "Method for cluster determination" },
1482 { "-minstruct", FALSE, etINT, {&minstruct},
1483 "Minimum number of structures in cluster for coloring in the [REF].xpm[ref] file" },
1484 { "-binary", FALSE, etBOOL, {&bBinary},
1485 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1486 "is given by [TT]-cutoff[tt]" },
1487 { "-M", FALSE, etINT, {&M},
1488 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1489 "0 is use cutoff" },
1490 { "-P", FALSE, etINT, {&P},
1491 "Number of identical nearest neighbors required to form a cluster" },
1492 { "-seed", FALSE, etINT, {&seed},
1493 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1494 { "-niter", FALSE, etINT, {&niter},
1495 "Number of iterations for MC" },
1496 { "-nrandom", FALSE, etINT, {&nrandom},
1497 "The first iterations for MC may be done complete random, to shuffle the frames" },
1498 { "-kT", FALSE, etREAL, {&kT},
1499 "Boltzmann weighting factor for Monte Carlo optimization "
1500 "(zero turns off uphill steps)" },
1501 { "-pbc", FALSE, etBOOL,
1502 { &bPBC }, "PBC check" }
1505 { efTRX, "-f", NULL, ffOPTRD },
1506 { efTPS, "-s", NULL, ffOPTRD },
1507 { efNDX, NULL, NULL, ffOPTRD },
1508 { efXPM, "-dm", "rmsd", ffOPTRD },
1509 { efXPM, "-om", "rmsd-raw", ffWRITE },
1510 { efXPM, "-o", "rmsd-clust", ffWRITE },
1511 { efLOG, "-g", "cluster", ffWRITE },
1512 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1513 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1514 { efXVG, "-conv", "mc-conv", ffOPTWR },
1515 { efXVG, "-sz", "clust-size", ffOPTWR},
1516 { efXPM, "-tr", "clust-trans", ffOPTWR},
1517 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1518 { efXVG, "-clid", "clust-id", ffOPTWR},
1519 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1521 #define NFILE asize(fnm)
1523 if (!parse_common_args(&argc, argv,
1524 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
1525 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1532 bReadMat = opt2bSet("-dm", NFILE, fnm);
1533 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1534 if (opt2parg_bSet("-av", asize(pa), pa) ||
1535 opt2parg_bSet("-wcl", asize(pa), pa) ||
1536 opt2parg_bSet("-nst", asize(pa), pa) ||
1537 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1538 opt2bSet("-cl", NFILE, fnm) )
1540 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1546 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1549 "\nWarning: assuming the time unit in %s is %s\n",
1550 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1552 if (trx_out_fn && !bReadTraj)
1554 fprintf(stderr, "\nWarning: "
1555 "cannot write cluster structures without reading trajectory\n"
1556 " ignoring option -cl %s\n", trx_out_fn);
1560 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1566 gmx_fatal(FARGS, "Invalid method");
1569 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1570 method == m_gromos );
1573 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1575 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1576 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1578 /* check input and write parameters to log file */
1579 bUseRmsdCut = FALSE;
1580 if (method == m_jarvis_patrick)
1582 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1583 if ((M < 0) || (M == 1))
1585 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1589 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1596 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1600 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1605 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1608 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1610 else /* method != m_jarvis */
1612 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1614 if (bUseRmsdCut && method != m_jarvis_patrick)
1616 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1618 if (method == m_monte_carlo)
1620 fprintf(log, "Using %d iterations\n", niter);
1625 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1631 /* don't read mass-database as masses (and top) are not used */
1632 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1636 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1639 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1640 bReadMat ? "" : " and RMSD calculation");
1641 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1642 1, &ifsize, &fitidx, &grpname);
1645 fprintf(stderr, "\nSelect group for output:\n");
1646 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1647 1, &iosize, &outidx, &grpname);
1648 /* merge and convert both index groups: */
1649 /* first copy outidx to index. let outidx refer to elements in index */
1650 snew(index, iosize);
1652 for (i = 0; i < iosize; i++)
1654 index[i] = outidx[i];
1657 /* now lookup elements from fitidx in index, add them if necessary
1658 and also let fitidx refer to elements in index */
1659 for (i = 0; i < ifsize; i++)
1662 while (j < isize && index[j] != fitidx[i])
1668 /* slow this way, but doesn't matter much */
1670 srenew(index, isize);
1672 index[j] = fitidx[i];
1676 else /* !trx_out_fn */
1680 for (i = 0; i < ifsize; i++)
1682 index[i] = fitidx[i];
1690 /* Loop over first coordinate file */
1691 fn = opt2fn("-f", NFILE, fnm);
1693 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1694 output_env_conv_times(oenv, nf, time);
1695 if (!bRMSdist || bAnalyze)
1697 /* Center all frames on zero */
1699 for (i = 0; i < ifsize; i++)
1701 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1705 for (i = 0; i < nf; i++)
1707 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1713 gmx_rmpbc_done(gpbc);
1719 fprintf(stderr, "Reading rms distance matrix ");
1720 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1721 fprintf(stderr, "\n");
1722 if (readmat[0].nx != readmat[0].ny)
1724 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1725 readmat[0].nx, readmat[0].ny);
1727 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1729 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1730 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1735 time = readmat[0].axis_x;
1736 time_invfac = output_env_get_time_invfactor(oenv);
1737 for (i = 0; i < nf; i++)
1739 time[i] *= time_invfac;
1742 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1743 convert_mat(&(readmat[0]), rms);
1745 nlevels = readmat[0].nmap;
1747 else /* !bReadMat */
1749 rms = init_mat(nf, method == m_diagonalize);
1750 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1753 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1754 /* Initialize work array */
1756 for (i1 = 0; i1 < nf; i1++)
1758 for (i2 = i1+1; i2 < nf; i2++)
1760 for (i = 0; i < isize; i++)
1762 copy_rvec(xx[i1][i], x1[i]);
1766 do_fit(isize, mass, xx[i2], x1);
1768 rmsd = rmsdev(isize, mass, xx[i2], x1);
1769 set_mat_entry(rms, i1, i2, rmsd);
1771 nrms -= (gmx_int64_t) (nf-i1-1);
1772 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1778 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1780 /* Initiate work arrays */
1783 for (i = 0; (i < isize); i++)
1788 for (i1 = 0; i1 < nf; i1++)
1790 calc_dist(isize, xx[i1], d1);
1791 for (i2 = i1+1; (i2 < nf); i2++)
1793 calc_dist(isize, xx[i2], d2);
1794 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1797 fprintf(stderr, "\r# RMSD calculations left: " "%"GMX_PRId64 " ", nrms);
1799 /* Clean up work arrays */
1800 for (i = 0; (i < isize); i++)
1808 fprintf(stderr, "\n\n");
1810 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1811 rms->minrms, rms->maxrms);
1812 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1813 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1814 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1815 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1817 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1818 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1820 if (bAnalyze && (rmsmin < rms->minrms) )
1822 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1823 rmsmin, rms->minrms);
1825 if (bAnalyze && (rmsmin > rmsdcut) )
1827 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1831 /* Plot the rmsd distribution */
1832 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1836 for (i1 = 0; (i1 < nf); i1++)
1838 for (i2 = 0; (i2 < nf); i2++)
1840 if (rms->mat[i1][i2] < rmsdcut)
1842 rms->mat[i1][i2] = 0;
1846 rms->mat[i1][i2] = 1;
1856 /* Now sort the matrix and write it out again */
1857 gather(rms, rmsdcut, &clust);
1860 /* Do a diagonalization */
1861 snew(eigenvalues, nf);
1862 snew(eigenvectors, nf*nf);
1863 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1864 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1865 sfree(eigenvectors);
1867 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1868 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1869 for (i = 0; (i < nf); i++)
1871 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1876 orig = init_mat(rms->nn, FALSE);
1878 copy_t_mat(orig, rms);
1879 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1880 opt2fn_null("-conv", NFILE, fnm), oenv);
1882 case m_jarvis_patrick:
1883 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1886 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1889 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1892 if (method == m_monte_carlo || method == m_diagonalize)
1894 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1902 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1906 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1908 init_t_atoms(&useatoms, isize, FALSE);
1909 snew(usextps, isize);
1910 useatoms.resinfo = top.atoms.resinfo;
1911 for (i = 0; i < isize; i++)
1913 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1914 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1915 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1);
1916 copy_rvec(xtps[index[i]], usextps[i]);
1918 useatoms.nr = isize;
1919 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1920 ifsize, fitidx, iosize, outidx,
1921 bReadTraj ? trx_out_fn : NULL,
1922 opt2fn_null("-sz", NFILE, fnm),
1923 opt2fn_null("-tr", NFILE, fnm),
1924 opt2fn_null("-ntr", NFILE, fnm),
1925 opt2fn_null("-clid", NFILE, fnm),
1926 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1927 rlo_bot, rhi_bot, oenv);
1931 if (bBinary && !bAnalyze)
1933 /* Make the clustering visible */
1934 for (i2 = 0; (i2 < nf); i2++)
1936 for (i1 = i2+1; (i1 < nf); i1++)
1938 if (rms->mat[i1][i2])
1940 rms->mat[i1][i2] = rms->maxrms;
1946 fp = opt2FILE("-o", NFILE, fnm, "w");
1947 fprintf(stderr, "Writing rms distance/clustering matrix ");
1950 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1951 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1952 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1956 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1957 sprintf(title, "RMS%sDeviation / Cluster Index",
1958 bRMSdist ? " Distance " : " ");
1961 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1962 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1963 rlo_top, rhi_top, 0.0, (real) ncluster,
1964 &ncluster, TRUE, rlo_bot, rhi_bot);
1968 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1969 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1970 rlo_top, rhi_top, &nlevels);
1973 fprintf(stderr, "\n");
1977 fp = opt2FILE("-om", NFILE, fnm, "w");
1978 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1979 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1980 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1981 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1982 rlo_top, rhi_top, &nlevels);
1987 /* now show what we've done */
1988 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1989 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1990 if (method == m_diagonalize)
1992 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1994 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1997 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1998 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1999 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
2001 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);