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.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/matio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/cmat.h"
51 #include "gromacs/gmxana/gmx_ana.h"
52 #include "gromacs/legacyheaders/macros.h"
53 #include "gromacs/legacyheaders/typedefs.h"
54 #include "gromacs/legacyheaders/viewit.h"
55 #include "gromacs/linearalgebra/eigensolver.h"
56 #include "gromacs/math/do_fit.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/pbcutil/pbc.h"
59 #include "gromacs/pbcutil/rmpbc.h"
60 #include "gromacs/random/random.h"
61 #include "gromacs/topology/index.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/smalloc.h"
67 /* print to two file pointers at once (i.e. stderr and log) */
69 void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
71 fprintf(fp1, "%s", buf);
72 fprintf(fp2, "%s", buf);
75 /* just print a prepared buffer to fp1 and fp2 */
77 void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
79 lo_ffprintf(fp1, fp2, buf);
82 /* prepare buffer with one argument, then print to fp1 and fp2 */
84 void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
86 sprintf(buf, fmt, arg);
87 lo_ffprintf(fp1, fp2, buf);
90 /* prepare buffer with one argument, then print to fp1 and fp2 */
92 void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
94 sprintf(buf, fmt, arg);
95 lo_ffprintf(fp1, fp2, buf);
98 /* prepare buffer with one argument, then print to fp1 and fp2 */
100 void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
102 sprintf(buf, fmt, arg);
103 lo_ffprintf(fp1, fp2, buf);
106 /* prepare buffer with two arguments, then print to fp1 and fp2 */
108 void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
110 sprintf(buf, fmt, arg1, arg2);
111 lo_ffprintf(fp1, fp2, buf);
114 /* prepare buffer with two arguments, then print to fp1 and fp2 */
116 void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
118 sprintf(buf, fmt, arg1, arg2);
119 lo_ffprintf(fp1, fp2, buf);
122 /* prepare buffer with two arguments, then print to fp1 and fp2 */
124 void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
126 sprintf(buf, fmt, arg1, arg2);
127 lo_ffprintf(fp1, fp2, buf);
140 void cp_index(int nn, int from[], int to[])
144 for (i = 0; (i < nn); i++)
150 void mc_optimize(FILE *log, t_mat *m, real *time,
151 int maxiter, int nrandom,
153 const char *conv, output_env_t oenv)
156 real ecur, enext, emin, prob, enorm;
157 int i, j, iswap, jswap, nn, nuphill = 0;
163 fprintf(stderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
166 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
167 printf("by reordering the frames to minimize the path between the two structures\n");
168 printf("that have the largest pairwise RMSD.\n");
171 enorm = m->mat[0][0];
172 for (i = 0; (i < m->n1); i++)
174 for (j = 0; (j < m->nn); j++)
176 if (m->mat[i][j] > enorm)
178 enorm = m->mat[i][j];
184 if ((iswap == -1) || (jswap == -1))
186 fprintf(stderr, "Matrix contains identical values in all fields\n");
189 swap_rows(m, 0, iswap);
190 swap_rows(m, m->n1-1, jswap);
191 emin = ecur = mat_energy(m);
192 printf("Largest distance %g between %d and %d. Energy: %g.\n",
193 enorm, iswap, jswap, emin);
195 rng = gmx_rng_init(seed);
198 /* Initiate and store global minimum */
199 minimum = init_mat(nn, m->b1D);
201 copy_t_mat(minimum, m);
205 fp = xvgropen(conv, "Convergence of the MC optimization",
206 "Energy", "Step", oenv);
208 for (i = 0; (i < maxiter); i++)
210 /* Generate new swapping candidates */
213 iswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
214 jswap = static_cast<int>(1+(nn-2)*gmx_rng_uniform_real(rng));
216 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
218 /* Apply swap and compute energy */
219 swap_rows(m, iswap, jswap);
220 enext = mat_energy(m);
222 /* Compute probability */
224 if ((enext < ecur) || (i < nrandom))
229 /* Store global minimum */
230 copy_t_mat(minimum, m);
236 /* Try Monte Carlo step */
237 prob = std::exp(-(enext-ecur)/(enorm*kT));
240 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
247 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
248 i, iswap, jswap, enext, prob);
251 fprintf(fp, "%6d %10g\n", i, enext);
257 swap_rows(m, jswap, iswap);
260 fprintf(log, "%d uphill steps were taken during optimization\n",
263 /* Now swap the matrix to get it into global minimum mode */
264 copy_t_mat(m, minimum);
266 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
267 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
268 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
269 for (i = 0; (i < m->nn); i++)
271 fprintf(log, "%10g %5d %10g\n",
274 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
283 static void calc_dist(int nind, rvec x[], real **d)
289 for (i = 0; (i < nind-1); i++)
292 for (j = i+1; (j < nind); j++)
294 /* Should use pbc_dx when analysing multiple molecueles,
295 * but the box is not stored for every frame.
297 rvec_sub(xi, x[j], dx);
303 static real rms_dist(int isize, real **d, real **d_r)
309 for (i = 0; (i < isize-1); i++)
311 for (j = i+1; (j < isize); j++)
313 r = d[i][j]-d_r[i][j];
317 r2 /= (isize*(isize-1))/2;
319 return std::sqrt(r2);
322 static int rms_dist_comp(const void *a, const void *b)
329 if (da->dist - db->dist < 0)
333 else if (da->dist - db->dist > 0)
340 static int clust_id_comp(const void *a, const void *b)
347 return da->clust - db->clust;
350 static int nrnb_comp(const void *a, const void *b)
357 /* return the b-a, we want highest first */
358 return db->nr - da->nr;
361 void gather(t_mat *m, real cutoff, t_clusters *clust)
365 int i, j, k, nn, cid, n1, diff;
368 /* First we sort the entries in the RMSD matrix */
372 for (i = k = 0; (i < n1); i++)
374 for (j = i+1; (j < n1); j++, k++)
378 d[k].dist = m->mat[i][j];
383 gmx_incons("gather algortihm");
385 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
387 /* Now we make a cluster index for all of the conformations */
390 /* Now we check the closest structures, and equalize their cluster numbers */
391 fprintf(stderr, "Linking structures ");
394 fprintf(stderr, "*");
396 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
398 diff = c[d[k].j].clust - c[d[k].i].clust;
404 c[d[k].j].clust = c[d[k].i].clust;
408 c[d[k].i].clust = c[d[k].j].clust;
414 fprintf(stderr, "\nSorting and renumbering clusters\n");
415 /* Sort on cluster number */
416 qsort(c, n1, sizeof(c[0]), clust_id_comp);
418 /* Renumber clusters */
420 for (k = 1; k < n1; k++)
422 if (c[k].clust != c[k-1].clust)
435 for (k = 0; (k < n1); k++)
437 fprintf(debug, "Cluster index for conformation %d: %d\n",
438 c[k].conf, c[k].clust);
442 for (k = 0; k < n1; k++)
444 clust->cl[c[k].conf] = c[k].clust;
451 gmx_bool jp_same(int **nnb, int i, int j, int P)
457 for (k = 0; nnb[i][k] >= 0; k++)
459 bIn = bIn || (nnb[i][k] == j);
467 for (k = 0; nnb[j][k] >= 0; k++)
469 bIn = bIn || (nnb[j][k] == i);
477 for (ii = 0; nnb[i][ii] >= 0; ii++)
479 for (jj = 0; nnb[j][jj] >= 0; jj++)
481 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
491 static void jarvis_patrick(int n1, real **mat, int M, int P,
492 real rmsdcut, t_clusters *clust)
497 int i, j, k, cid, diff, maxval;
506 /* First we sort the entries in the RMSD matrix row by row.
507 * This gives us the nearest neighbor list.
511 for (i = 0; (i < n1); i++)
513 for (j = 0; (j < n1); j++)
516 row[j].dist = mat[i][j];
518 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
521 /* Put the M nearest neighbors in the list */
523 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
527 nnb[i][k] = row[j].j;
535 /* Put all neighbors nearer than rmsdcut in the list */
538 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
545 srenew(nnb[i], maxval);
547 nnb[i][k] = row[j].j;
553 srenew(nnb[i], maxval+1);
561 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
562 for (i = 0; (i < n1); i++)
564 fprintf(debug, "i:%5d nbs:", i);
565 for (j = 0; nnb[i][j] >= 0; j++)
567 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
569 fprintf(debug, "\n");
574 fprintf(stderr, "Linking structures ");
575 /* Use mcpy for temporary storage of booleans */
576 mcpy = mk_matrix(n1, n1, FALSE);
577 for (i = 0; i < n1; i++)
579 for (j = i+1; j < n1; j++)
581 mcpy[i][j] = jp_same(nnb, i, j, P);
586 fprintf(stderr, "*");
588 for (i = 0; i < n1; i++)
590 for (j = i+1; j < n1; j++)
594 diff = c[j].clust - c[i].clust;
600 c[j].clust = c[i].clust;
604 c[i].clust = c[j].clust;
613 fprintf(stderr, "\nSorting and renumbering clusters\n");
614 /* Sort on cluster number */
615 qsort(c, n1, sizeof(c[0]), clust_id_comp);
617 /* Renumber clusters */
619 for (k = 1; k < n1; k++)
621 if (c[k].clust != c[k-1].clust)
633 for (k = 0; k < n1; k++)
635 clust->cl[c[k].conf] = c[k].clust;
639 for (k = 0; (k < n1); k++)
641 fprintf(debug, "Cluster index for conformation %d: %d\n",
642 c[k].conf, c[k].clust);
646 done_matrix(n1, &mcpy);
649 for (i = 0; (i < n1); i++)
656 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
660 /* dump neighbor list */
661 fprintf(fp, "%s", title);
662 for (i = 0; (i < n1); i++)
664 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
665 for (j = 0; j < nnb[i].nr; j++)
667 fprintf(fp, "%5d", nnb[i].nb[j]);
673 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
677 int i, j, k, j1, maxval;
679 /* Put all neighbors nearer than rmsdcut in the list */
680 fprintf(stderr, "Making list of neighbors within cutoff ");
683 for (i = 0; (i < n1); i++)
687 /* put all neighbors within cut-off in list */
688 for (j = 0; j < n1; j++)
690 if (mat[i][j] < rmsdcut)
695 srenew(nnb[i].nb, maxval);
701 /* store nr of neighbors, we'll need that */
703 if (i%(1+n1/100) == 0)
705 fprintf(stderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
708 fprintf(stderr, "%3d%%\n", 100);
711 /* sort neighbor list on number of neighbors, largest first */
712 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
716 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
719 /* turn first structure with all its neighbors (largest) into cluster
720 remove them from pool of structures and repeat for all remaining */
721 fprintf(stderr, "Finding clusters %4d", 0);
722 /* cluster id's start at 1: */
726 /* set cluster id (k) for first item in neighborlist */
727 for (j = 0; j < nnb[0].nr; j++)
729 clust->cl[nnb[0].nb[j]] = k;
735 /* adjust number of neighbors for others, taking removals into account: */
736 for (i = 1; i < n1 && nnb[i].nr; i++)
739 for (j = 0; j < nnb[i].nr; j++)
741 /* if this neighbor wasn't removed */
742 if (clust->cl[nnb[i].nb[j]] == 0)
744 /* shift the rest (j1<=j) */
745 nnb[i].nb[j1] = nnb[i].nb[j];
750 /* now j1 is the new number of neighbors */
753 /* sort again on nnb[].nr, because we have new # neighbors: */
754 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
755 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
757 fprintf(stderr, "\b\b\b\b%4d", k);
761 fprintf(stderr, "\n");
765 fprintf(debug, "Clusters (%d):\n", k);
766 for (i = 0; i < n1; i++)
768 fprintf(debug, " %3d", clust->cl[i]);
770 fprintf(debug, "\n");
776 rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
777 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
782 int i, i0, j, max_nf;
790 natom = read_first_x(oenv, &status, fn, &t, &x, box);
797 gmx_rmpbc(gpbc, natom, box, x);
803 srenew(*time, max_nf);
808 /* Store only the interesting atoms */
809 for (j = 0; (j < isize); j++)
811 copy_rvec(x[index[j]], xx[i0][j]);
818 while (read_next_x(oenv, status, &t, x, box));
819 fprintf(stderr, "Allocated %lu bytes for frames\n",
820 (unsigned long) (max_nf*isize*sizeof(**xx)));
821 fprintf(stderr, "Read %d frames from trajectory %s\n", i0, fn);
828 static int plot_clusters(int nf, real **mat, t_clusters *clust,
831 int i, j, ncluster, ci;
832 int *cl_id, *nstruct, *strind;
837 for (i = 0; i < nf; i++)
840 cl_id[i] = clust->cl[i];
844 for (i = 0; i < nf; i++)
846 if (nstruct[i] >= minstruct)
849 for (j = 0; (j < nf); j++)
853 strind[j] = ncluster;
859 fprintf(stderr, "There are %d clusters with at least %d conformations\n",
860 ncluster, minstruct);
862 for (i = 0; (i < nf); i++)
865 for (j = 0; j < i; j++)
867 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
869 /* color different clusters with different colors, as long as
870 we don't run out of colors */
871 mat[i][j] = strind[i];
886 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
890 for (i = 0; i < nf; i++)
892 for (j = 0; j < i; j++)
894 if (clust->cl[i] == clust->cl[j])
906 static char *parse_filename(const char *fn, int maxnr)
913 if (std::strchr(fn, '%'))
915 gmx_fatal(FARGS, "will not number filename %s containing '%c'", fn, '%');
917 /* number of digits needed in numbering */
918 i = static_cast<int>((std::log(static_cast<real>(maxnr))/std::log(10.0)) + 1);
919 /* split fn and ext */
920 ext = std::strrchr(fn, '.');
923 gmx_fatal(FARGS, "cannot separate extension in filename %s", fn);
926 /* insert e.g. '%03d' between fn and ext */
927 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
928 snew(fnout, std::strlen(buf)+1);
929 std::strcpy(fnout, buf);
934 static void ana_trans(t_clusters *clust, int nf,
935 const char *transfn, const char *ntransfn, FILE *log,
936 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
941 int i, ntranst, maxtrans;
944 snew(ntrans, clust->ncl);
945 snew(trans, clust->ncl);
946 snew(axis, clust->ncl);
947 for (i = 0; i < clust->ncl; i++)
950 snew(trans[i], clust->ncl);
954 for (i = 1; i < nf; i++)
956 if (clust->cl[i] != clust->cl[i-1])
959 ntrans[clust->cl[i-1]-1]++;
960 ntrans[clust->cl[i]-1]++;
961 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
962 maxtrans = static_cast<int>(std::max(static_cast<real>(maxtrans), trans[clust->cl[i]-1][clust->cl[i-1]-1]));
965 ffprintf_dd(stderr, log, buf, "Counted %d transitions in total, "
966 "max %d between two specific clusters\n", ntranst, maxtrans);
969 fp = gmx_ffopen(transfn, "w");
970 i = std::min(maxtrans+1, 80);
971 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
972 "from cluster", "to cluster",
973 clust->ncl, clust->ncl, axis, axis, trans,
974 0, maxtrans, rlo, rhi, &i);
979 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
981 for (i = 0; i < clust->ncl; i++)
983 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
988 for (i = 0; i < clust->ncl; i++)
996 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
997 int natom, t_atoms *atoms, rvec *xtps,
998 real *mass, rvec **xx, real *time,
999 int ifsize, atom_id *fitidx,
1000 int iosize, atom_id *outidx,
1001 const char *trxfn, const char *sizefn,
1002 const char *transfn, const char *ntransfn,
1003 const char *clustidfn, gmx_bool bAverage,
1004 int write_ncl, int write_nst, real rmsmin,
1005 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1006 const output_env_t oenv)
1008 FILE *size_fp = NULL;
1009 char buf[STRLEN], buf1[40], buf2[40], buf3[40], *trxsfn;
1010 t_trxstatus *trxout = NULL;
1011 t_trxstatus *trxsout = NULL;
1012 int i, i1, cl, nstr, *structure, first = 0, midstr;
1013 gmx_bool *bWrite = NULL;
1014 real r, clrmsd, midrmsd;
1020 ffprintf_d(stderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1024 /* do we write all structures? */
1027 trxsfn = parse_filename(trxfn, std::max(write_ncl, clust->ncl));
1030 ffprintf_ss(stderr, log, buf, "Writing %s structure for each cluster to %s\n",
1031 bAverage ? "average" : "middle", trxfn);
1034 /* find out what we want to tell the user:
1035 Writing [all structures|structures with rmsd > %g] for
1036 {all|first %d} clusters {with more than %d structures} to %s */
1039 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1043 sprintf(buf1, "all structures");
1045 buf2[0] = buf3[0] = '\0';
1046 if (write_ncl >= clust->ncl)
1050 sprintf(buf2, "all ");
1055 sprintf(buf2, "the first %d ", write_ncl);
1059 sprintf(buf3, " with more than %d structures", write_nst);
1061 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1062 ffprintf(stderr, log, buf);
1065 /* Prepare a reference structure for the orientation of the clusters */
1068 reset_x(ifsize, fitidx, natom, NULL, xtps, mass);
1070 trxout = open_trx(trxfn, "w");
1071 /* Calculate the average structure in each cluster, *
1072 * all structures are fitted to the first struture of the cluster */
1076 if (transfn || ntransfn)
1078 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1083 FILE *fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1084 if (output_env_get_print_xvgr_codes(oenv))
1086 fprintf(fp, "@ s0 symbol 2\n");
1087 fprintf(fp, "@ s0 symbol size 0.2\n");
1088 fprintf(fp, "@ s0 linestyle 0\n");
1090 for (i = 0; i < nf; i++)
1092 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1098 size_fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1099 if (output_env_get_print_xvgr_codes(oenv))
1101 fprintf(size_fp, "@g%d type %s\n", 0, "bar");
1104 snew(structure, nf);
1105 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1106 "cl.", "#st", "rmsd", "middle", "rmsd");
1107 for (cl = 1; cl <= clust->ncl; cl++)
1109 /* prepare structures (fit, middle, average) */
1112 for (i = 0; i < natom; i++)
1118 for (i1 = 0; i1 < nf; i1++)
1120 if (clust->cl[i1] == cl)
1122 structure[nstr] = i1;
1124 if (trxfn && (bAverage || write_ncl) )
1128 reset_x(ifsize, fitidx, natom, NULL, xx[i1], mass);
1136 do_fit(natom, mass, xx[first], xx[i1]);
1140 for (i = 0; i < natom; i++)
1142 rvec_inc(xav[i], xx[i1][i]);
1150 fprintf(size_fp, "%8d %8d\n", cl, nstr);
1155 for (i1 = 0; i1 < nstr; i1++)
1160 for (i = 0; i < nstr; i++)
1164 r += rmsd[structure[i]][structure[i1]];
1168 r += rmsd[structure[i1]][structure[i]];
1175 midstr = structure[i1];
1182 /* dump cluster info to logfile */
1185 sprintf(buf1, "%6.3f", clrmsd);
1190 sprintf(buf2, "%5.3f", midrmsd);
1198 sprintf(buf1, "%5s", "");
1199 sprintf(buf2, "%5s", "");
1201 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1202 for (i = 0; i < nstr; i++)
1204 if ((i % 7 == 0) && i)
1206 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1213 fprintf(log, "%s %6g", buf, time[i1]);
1217 /* write structures to trajectory file(s) */
1222 for (i = 0; i < nstr; i++)
1227 if (cl < write_ncl+1 && nstr > write_nst)
1229 /* Dump all structures for this cluster */
1230 /* generate numbered filename (there is a %d in trxfn!) */
1231 sprintf(buf, trxsfn, cl);
1232 trxsout = open_trx(buf, "w");
1233 for (i = 0; i < nstr; i++)
1238 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1242 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1248 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1249 xx[structure[i]], NULL, NULL);
1254 /* Dump the average structure for this cluster */
1257 for (i = 0; i < natom; i++)
1259 svmul(1.0/nstr, xav[i], xav[i]);
1264 for (i = 0; i < natom; i++)
1266 copy_rvec(xx[midstr][i], xav[i]);
1270 reset_x(ifsize, fitidx, natom, NULL, xav, mass);
1275 do_fit(natom, mass, xtps, xav);
1277 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL, NULL);
1302 static void convert_mat(t_matrix *mat, t_mat *rms)
1307 matrix2real(mat, rms->mat);
1308 /* free input xpm matrix data */
1309 for (i = 0; i < mat->nx; i++)
1311 sfree(mat->matrix[i]);
1315 for (i = 0; i < mat->nx; i++)
1317 for (j = i; j < mat->nx; j++)
1319 rms->sumrms += rms->mat[i][j];
1320 rms->maxrms = std::max(rms->maxrms, rms->mat[i][j]);
1323 rms->minrms = std::min(rms->minrms, rms->mat[i][j]);
1330 int gmx_cluster(int argc, char *argv[])
1332 const char *desc[] = {
1333 "[THISMODULE] can cluster structures using several different methods.",
1334 "Distances between structures can be determined from a trajectory",
1335 "or read from an [REF].xpm[ref] matrix file with the [TT]-dm[tt] option.",
1336 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1337 "can be used to define the distance between structures.[PAR]",
1339 "single linkage: add a structure to a cluster when its distance to any",
1340 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1342 "Jarvis Patrick: add a structure to a cluster when this structure",
1343 "and a structure in the cluster have each other as neighbors and",
1344 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1345 "of a structure are the M closest structures or all structures within",
1346 "[TT]cutoff[tt].[PAR]",
1348 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1349 "the order of the frames is using the smallest possible increments.",
1350 "With this it is possible to make a smooth animation going from one",
1351 "structure to another with the largest possible (e.g.) RMSD between",
1352 "them, however the intermediate steps should be as small as possible.",
1353 "Applications could be to visualize a potential of mean force",
1354 "ensemble of simulations or a pulling simulation. Obviously the user",
1355 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1356 "The final result can be inspect visually by looking at the matrix",
1357 "[REF].xpm[ref] file, which should vary smoothly from bottom to top.[PAR]",
1359 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1361 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1362 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1363 "Count number of neighbors using cut-off, take structure with",
1364 "largest number of neighbors with all its neighbors as cluster",
1365 "and eliminate it from the pool of clusters. Repeat for remaining",
1366 "structures in pool.[PAR]",
1368 "When the clustering algorithm assigns each structure to exactly one",
1369 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1370 "file is supplied, the structure with",
1371 "the smallest average distance to the others or the average structure",
1372 "or all structures for each cluster will be written to a trajectory",
1373 "file. When writing all structures, separate numbered files are made",
1374 "for each cluster.[PAR]",
1376 "Two output files are always written:",
1378 " * [TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1379 " and a graphical depiction of the clusters in the lower right half",
1380 " When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1381 " when two structures are in the same cluster.",
1382 " When [TT]-minstruct[tt] > 1 different colors will be used for each",
1384 " * [TT]-g[tt] writes information on the options used and a detailed list",
1385 " of all clusters and their members.",
1388 "Additionally, a number of optional output files can be written:",
1390 " * [TT]-dist[tt] writes the RMSD distribution.",
1391 " * [TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1392 " diagonalization.",
1393 " * [TT]-sz[tt] writes the cluster sizes.",
1394 " * [TT]-tr[tt] writes a matrix of the number transitions between",
1396 " * [TT]-ntr[tt] writes the total number of transitions to or from",
1398 " * [TT]-clid[tt] writes the cluster number as a function of time.",
1399 " * [TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1400 " structure of each cluster or writes numbered files with cluster members",
1401 " for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1402 " [TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1403 " structure with the smallest average RMSD from all other structures",
1408 int nf, i, i1, i2, j;
1409 gmx_int64_t nrms = 0;
1412 rvec *xtps, *usextps, *x1, **xx = NULL;
1413 const char *fn, *trx_out_fn;
1415 t_mat *rms, *orig = NULL;
1420 t_matrix *readmat = NULL;
1423 int isize = 0, ifsize = 0, iosize = 0;
1424 atom_id *index = NULL, *fitidx = NULL, *outidx = NULL;
1426 real rmsd, **d1, **d2, *time = NULL, time_invfac, *mass = NULL;
1427 char buf[STRLEN], buf1[80], title[STRLEN];
1428 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE, bReadMat, bReadTraj, bPBC = TRUE;
1430 int method, ncluster = 0;
1431 static const char *methodname[] = {
1432 NULL, "linkage", "jarvis-patrick", "monte-carlo",
1433 "diagonalization", "gromos", NULL
1436 m_null, m_linkage, m_jarvis_patrick,
1437 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1439 /* Set colors for plotting: white = zero RMS, black = maximum */
1440 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1441 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1442 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1443 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1444 static int nlevels = 40, skip = 1;
1445 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1446 gmx_bool bRMSdist = FALSE, bBinary = FALSE, bAverage = FALSE, bFit = TRUE;
1447 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1448 static real kT = 1e-3;
1449 static int M = 10, P = 3;
1451 gmx_rmpbc_t gpbc = NULL;
1454 { "-dista", FALSE, etBOOL, {&bRMSdist},
1455 "Use RMSD of distances instead of RMS deviation" },
1456 { "-nlevels", FALSE, etINT, {&nlevels},
1457 "Discretize RMSD matrix in this number of levels" },
1458 { "-cutoff", FALSE, etREAL, {&rmsdcut},
1459 "RMSD cut-off (nm) for two structures to be neighbor" },
1460 { "-fit", FALSE, etBOOL, {&bFit},
1461 "Use least squares fitting before RMSD calculation" },
1462 { "-max", FALSE, etREAL, {&scalemax},
1463 "Maximum level in RMSD matrix" },
1464 { "-skip", FALSE, etINT, {&skip},
1465 "Only analyze every nr-th frame" },
1466 { "-av", FALSE, etBOOL, {&bAverage},
1467 "Write average iso middle structure for each cluster" },
1468 { "-wcl", FALSE, etINT, {&write_ncl},
1469 "Write the structures for this number of clusters to numbered files" },
1470 { "-nst", FALSE, etINT, {&write_nst},
1471 "Only write all structures if more than this number of structures per cluster" },
1472 { "-rmsmin", FALSE, etREAL, {&rmsmin},
1473 "minimum rms difference with rest of cluster for writing structures" },
1474 { "-method", FALSE, etENUM, {methodname},
1475 "Method for cluster determination" },
1476 { "-minstruct", FALSE, etINT, {&minstruct},
1477 "Minimum number of structures in cluster for coloring in the [REF].xpm[ref] file" },
1478 { "-binary", FALSE, etBOOL, {&bBinary},
1479 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1480 "is given by [TT]-cutoff[tt]" },
1481 { "-M", FALSE, etINT, {&M},
1482 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1483 "0 is use cutoff" },
1484 { "-P", FALSE, etINT, {&P},
1485 "Number of identical nearest neighbors required to form a cluster" },
1486 { "-seed", FALSE, etINT, {&seed},
1487 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1488 { "-niter", FALSE, etINT, {&niter},
1489 "Number of iterations for MC" },
1490 { "-nrandom", FALSE, etINT, {&nrandom},
1491 "The first iterations for MC may be done complete random, to shuffle the frames" },
1492 { "-kT", FALSE, etREAL, {&kT},
1493 "Boltzmann weighting factor for Monte Carlo optimization "
1494 "(zero turns off uphill steps)" },
1495 { "-pbc", FALSE, etBOOL,
1496 { &bPBC }, "PBC check" }
1499 { efTRX, "-f", NULL, ffOPTRD },
1500 { efTPS, "-s", NULL, ffOPTRD },
1501 { efNDX, NULL, NULL, ffOPTRD },
1502 { efXPM, "-dm", "rmsd", ffOPTRD },
1503 { efXPM, "-om", "rmsd-raw", ffWRITE },
1504 { efXPM, "-o", "rmsd-clust", ffWRITE },
1505 { efLOG, "-g", "cluster", ffWRITE },
1506 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
1507 { efXVG, "-ev", "rmsd-eig", ffOPTWR },
1508 { efXVG, "-conv", "mc-conv", ffOPTWR },
1509 { efXVG, "-sz", "clust-size", ffOPTWR},
1510 { efXPM, "-tr", "clust-trans", ffOPTWR},
1511 { efXVG, "-ntr", "clust-trans", ffOPTWR},
1512 { efXVG, "-clid", "clust-id", ffOPTWR},
1513 { efTRX, "-cl", "clusters.pdb", ffOPTWR }
1515 #define NFILE asize(fnm)
1517 if (!parse_common_args(&argc, argv,
1518 PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT,
1519 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
1526 bReadMat = opt2bSet("-dm", NFILE, fnm);
1527 bReadTraj = opt2bSet("-f", NFILE, fnm) || !bReadMat;
1528 if (opt2parg_bSet("-av", asize(pa), pa) ||
1529 opt2parg_bSet("-wcl", asize(pa), pa) ||
1530 opt2parg_bSet("-nst", asize(pa), pa) ||
1531 opt2parg_bSet("-rmsmin", asize(pa), pa) ||
1532 opt2bSet("-cl", NFILE, fnm) )
1534 trx_out_fn = opt2fn("-cl", NFILE, fnm);
1540 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1543 "\nWarning: assuming the time unit in %s is %s\n",
1544 opt2fn("-dm", NFILE, fnm), output_env_get_time_unit(oenv));
1546 if (trx_out_fn && !bReadTraj)
1548 fprintf(stderr, "\nWarning: "
1549 "cannot write cluster structures without reading trajectory\n"
1550 " ignoring option -cl %s\n", trx_out_fn);
1554 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1560 gmx_fatal(FARGS, "Invalid method");
1563 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1564 method == m_gromos );
1567 log = ftp2FILE(efLOG, NFILE, fnm, "w");
1569 fprintf(stderr, "Using %s method for clustering\n", methodname[0]);
1570 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1572 /* check input and write parameters to log file */
1573 bUseRmsdCut = FALSE;
1574 if (method == m_jarvis_patrick)
1576 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa), pa);
1577 if ((M < 0) || (M == 1))
1579 gmx_fatal(FARGS, "M (%d) must be 0 or larger than 1", M);
1583 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1590 gmx_fatal(FARGS, "Number of neighbors required (P) must be less than M");
1594 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1599 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1602 ffprintf_s(stderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1604 else /* method != m_jarvis */
1606 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1608 if (bUseRmsdCut && method != m_jarvis_patrick)
1610 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1612 if (method == m_monte_carlo)
1614 fprintf(log, "Using %d iterations\n", niter);
1619 gmx_fatal(FARGS, "skip (%d) should be >= 1", skip);
1625 /* don't read mass-database as masses (and top) are not used */
1626 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xtps, NULL, box,
1630 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1633 fprintf(stderr, "\nSelect group for least squares fit%s:\n",
1634 bReadMat ? "" : " and RMSD calculation");
1635 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1636 1, &ifsize, &fitidx, &grpname);
1639 fprintf(stderr, "\nSelect group for output:\n");
1640 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
1641 1, &iosize, &outidx, &grpname);
1642 /* merge and convert both index groups: */
1643 /* first copy outidx to index. let outidx refer to elements in index */
1644 snew(index, iosize);
1646 for (i = 0; i < iosize; i++)
1648 index[i] = outidx[i];
1651 /* now lookup elements from fitidx in index, add them if necessary
1652 and also let fitidx refer to elements in index */
1653 for (i = 0; i < ifsize; i++)
1656 while (j < isize && index[j] != fitidx[i])
1662 /* slow this way, but doesn't matter much */
1664 srenew(index, isize);
1666 index[j] = fitidx[i];
1670 else /* !trx_out_fn */
1674 for (i = 0; i < ifsize; i++)
1676 index[i] = fitidx[i];
1684 /* Loop over first coordinate file */
1685 fn = opt2fn("-f", NFILE, fnm);
1687 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1688 output_env_conv_times(oenv, nf, time);
1689 if (!bRMSdist || bAnalyze)
1691 /* Center all frames on zero */
1693 for (i = 0; i < ifsize; i++)
1695 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1699 for (i = 0; i < nf; i++)
1701 reset_x(ifsize, fitidx, isize, NULL, xx[i], mass);
1707 gmx_rmpbc_done(gpbc);
1713 fprintf(stderr, "Reading rms distance matrix ");
1714 read_xpm_matrix(opt2fn("-dm", NFILE, fnm), &readmat);
1715 fprintf(stderr, "\n");
1716 if (readmat[0].nx != readmat[0].ny)
1718 gmx_fatal(FARGS, "Matrix (%dx%d) is not square",
1719 readmat[0].nx, readmat[0].ny);
1721 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1723 gmx_fatal(FARGS, "Matrix size (%dx%d) does not match the number of "
1724 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1729 time = readmat[0].axis_x;
1730 time_invfac = output_env_get_time_invfactor(oenv);
1731 for (i = 0; i < nf; i++)
1733 time[i] *= time_invfac;
1736 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1737 convert_mat(&(readmat[0]), rms);
1739 nlevels = readmat[0].nmap;
1741 else /* !bReadMat */
1743 rms = init_mat(nf, method == m_diagonalize);
1744 nrms = (static_cast<gmx_int64_t>(nf)*static_cast<gmx_int64_t>(nf-1))/2;
1747 fprintf(stderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1748 /* Initialize work array */
1750 for (i1 = 0; i1 < nf; i1++)
1752 for (i2 = i1+1; i2 < nf; i2++)
1754 for (i = 0; i < isize; i++)
1756 copy_rvec(xx[i1][i], x1[i]);
1760 do_fit(isize, mass, xx[i2], x1);
1762 rmsd = rmsdev(isize, mass, xx[i2], x1);
1763 set_mat_entry(rms, i1, i2, rmsd);
1766 fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
1772 fprintf(stderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1774 /* Initiate work arrays */
1777 for (i = 0; (i < isize); i++)
1782 for (i1 = 0; i1 < nf; i1++)
1784 calc_dist(isize, xx[i1], d1);
1785 for (i2 = i1+1; (i2 < nf); i2++)
1787 calc_dist(isize, xx[i2], d2);
1788 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1791 fprintf(stderr, "\r# RMSD calculations left: " "%" GMX_PRId64 " ", nrms);
1793 /* Clean up work arrays */
1794 for (i = 0; (i < isize); i++)
1802 fprintf(stderr, "\n\n");
1804 ffprintf_gg(stderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1805 rms->minrms, rms->maxrms);
1806 ffprintf_g(stderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1807 ffprintf_d(stderr, log, buf, "Number of structures for matrix %d\n", nf);
1808 ffprintf_g(stderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1809 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1811 fprintf(stderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1812 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1814 if (bAnalyze && (rmsmin < rms->minrms) )
1816 fprintf(stderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1817 rmsmin, rms->minrms);
1819 if (bAnalyze && (rmsmin > rmsdcut) )
1821 fprintf(stderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1825 /* Plot the rmsd distribution */
1826 rmsd_distribution(opt2fn("-dist", NFILE, fnm), rms, oenv);
1830 for (i1 = 0; (i1 < nf); i1++)
1832 for (i2 = 0; (i2 < nf); i2++)
1834 if (rms->mat[i1][i2] < rmsdcut)
1836 rms->mat[i1][i2] = 0;
1840 rms->mat[i1][i2] = 1;
1850 /* Now sort the matrix and write it out again */
1851 gather(rms, rmsdcut, &clust);
1854 /* Do a diagonalization */
1855 snew(eigenvalues, nf);
1856 snew(eigenvectors, nf*nf);
1857 std::memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1858 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1859 sfree(eigenvectors);
1861 fp = xvgropen(opt2fn("-ev", NFILE, fnm), "RMSD matrix Eigenvalues",
1862 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1863 for (i = 0; (i < nf); i++)
1865 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1870 orig = init_mat(rms->nn, FALSE);
1872 copy_t_mat(orig, rms);
1873 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1874 opt2fn_null("-conv", NFILE, fnm), oenv);
1876 case m_jarvis_patrick:
1877 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1880 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1883 gmx_fatal(FARGS, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1886 if (method == m_monte_carlo || method == m_diagonalize)
1888 fprintf(stderr, "Energy of the matrix after clustering is %g.\n",
1896 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1900 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1902 init_t_atoms(&useatoms, isize, FALSE);
1903 snew(usextps, isize);
1904 useatoms.resinfo = top.atoms.resinfo;
1905 for (i = 0; i < isize; i++)
1907 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1908 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1909 useatoms.nres = std::max(useatoms.nres, useatoms.atom[i].resind+1);
1910 copy_rvec(xtps[index[i]], usextps[i]);
1912 useatoms.nr = isize;
1913 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1914 ifsize, fitidx, iosize, outidx,
1915 bReadTraj ? trx_out_fn : NULL,
1916 opt2fn_null("-sz", NFILE, fnm),
1917 opt2fn_null("-tr", NFILE, fnm),
1918 opt2fn_null("-ntr", NFILE, fnm),
1919 opt2fn_null("-clid", NFILE, fnm),
1920 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1921 rlo_bot, rhi_bot, oenv);
1925 if (bBinary && !bAnalyze)
1927 /* Make the clustering visible */
1928 for (i2 = 0; (i2 < nf); i2++)
1930 for (i1 = i2+1; (i1 < nf); i1++)
1932 if (rms->mat[i1][i2])
1934 rms->mat[i1][i2] = rms->maxrms;
1940 fp = opt2FILE("-o", NFILE, fnm, "w");
1941 fprintf(stderr, "Writing rms distance/clustering matrix ");
1944 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1945 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1946 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1950 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1951 sprintf(title, "RMS%sDeviation / Cluster Index",
1952 bRMSdist ? " Distance " : " ");
1955 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1956 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1957 rlo_top, rhi_top, 0.0, ncluster,
1958 &ncluster, TRUE, rlo_bot, rhi_bot);
1962 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1963 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1964 rlo_top, rhi_top, &nlevels);
1967 fprintf(stderr, "\n");
1971 fp = opt2FILE("-om", NFILE, fnm, "w");
1972 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1973 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1974 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1975 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1976 rlo_top, rhi_top, &nlevels);
1981 /* now show what we've done */
1982 do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
1983 do_view(oenv, opt2fn_null("-sz", NFILE, fnm), "-nxy");
1984 if (method == m_diagonalize)
1986 do_view(oenv, opt2fn_null("-ev", NFILE, fnm), "-nxy");
1988 do_view(oenv, opt2fn("-dist", NFILE, fnm), "-nxy");
1991 do_view(oenv, opt2fn_null("-tr", NFILE, fnm), "-nxy");
1992 do_view(oenv, opt2fn_null("-ntr", NFILE, fnm), "-nxy");
1993 do_view(oenv, opt2fn_null("-clid", NFILE, fnm), "-nxy");
1995 do_view(oenv, opt2fn_null("-conv", NFILE, fnm), NULL);