Bug Summary

File:gromacs/gmxana/gmx_cluster.c
Location:line 1288, column 13
Description:Value stored to 'r' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <math.h>
42#include <stdlib.h>
43#include <string.h>
44
45#include "macros.h"
46#include "gromacs/utility/smalloc.h"
47#include "typedefs.h"
48#include "gromacs/commandline/pargs.h"
49#include "gromacs/fileio/tpxio.h"
50#include "gromacs/fileio/trxio.h"
51#include "gromacs/utility/cstringutil.h"
52#include "gromacs/math/vec.h"
53#include "macros.h"
54#include "index.h"
55#include "gromacs/random/random.h"
56#include "pbc.h"
57#include "rmpbc.h"
58#include "gromacs/fileio/xvgr.h"
59#include "gromacs/utility/futil.h"
60#include "gromacs/fileio/matio.h"
61#include "cmat.h"
62#include "gromacs/fileio/trnio.h"
63#include "viewit.h"
64#include "gmx_ana.h"
65
66#include "gromacs/linearalgebra/eigensolver.h"
67#include "gromacs/math/do_fit.h"
68#include "gromacs/utility/fatalerror.h"
69
70/* print to two file pointers at once (i.e. stderr and log) */
71static gmx_inlineinline
72void lo_ffprintf(FILE *fp1, FILE *fp2, const char *buf)
73{
74 fprintf(fp1, "%s", buf);
75 fprintf(fp2, "%s", buf);
76}
77
78/* just print a prepared buffer to fp1 and fp2 */
79static gmx_inlineinline
80void ffprintf(FILE *fp1, FILE *fp2, const char *buf)
81{
82 lo_ffprintf(fp1, fp2, buf);
83}
84
85/* prepare buffer with one argument, then print to fp1 and fp2 */
86static gmx_inlineinline
87void ffprintf_d(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg)
88{
89 sprintf(buf, fmt, arg);
90 lo_ffprintf(fp1, fp2, buf);
91}
92
93/* prepare buffer with one argument, then print to fp1 and fp2 */
94static gmx_inlineinline
95void ffprintf_g(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg)
96{
97 sprintf(buf, fmt, arg);
98 lo_ffprintf(fp1, fp2, buf);
99}
100
101/* prepare buffer with one argument, then print to fp1 and fp2 */
102static gmx_inlineinline
103void ffprintf_s(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg)
104{
105 sprintf(buf, fmt, arg);
106 lo_ffprintf(fp1, fp2, buf);
107}
108
109/* prepare buffer with two arguments, then print to fp1 and fp2 */
110static gmx_inlineinline
111void ffprintf_dd(FILE *fp1, FILE *fp2, char *buf, const char *fmt, int arg1, int arg2)
112{
113 sprintf(buf, fmt, arg1, arg2);
114 lo_ffprintf(fp1, fp2, buf);
115}
116
117/* prepare buffer with two arguments, then print to fp1 and fp2 */
118static gmx_inlineinline
119void ffprintf_gg(FILE *fp1, FILE *fp2, char *buf, const char *fmt, real arg1, real arg2)
120{
121 sprintf(buf, fmt, arg1, arg2);
122 lo_ffprintf(fp1, fp2, buf);
123}
124
125/* prepare buffer with two arguments, then print to fp1 and fp2 */
126static gmx_inlineinline
127void ffprintf_ss(FILE *fp1, FILE *fp2, char *buf, const char *fmt, const char *arg1, const char *arg2)
128{
129 sprintf(buf, fmt, arg1, arg2);
130 lo_ffprintf(fp1, fp2, buf);
131}
132
133typedef struct {
134 int ncl;
135 int *cl;
136} t_clusters;
137
138typedef struct {
139 int nr;
140 int *nb;
141} t_nnb;
142
143void pr_energy(FILE *fp, real e)
144{
145 fprintf(fp, "Energy: %8.4f\n", e);
146}
147
148void cp_index(int nn, int from[], int to[])
149{
150 int i;
151
152 for (i = 0; (i < nn); i++)
153 {
154 to[i] = from[i];
155 }
156}
157
158void mc_optimize(FILE *log, t_mat *m, real *time,
159 int maxiter, int nrandom,
160 int seed, real kT,
161 const char *conv, output_env_t oenv)
162{
163 FILE *fp = NULL((void*)0);
164 real ecur, enext, emin, prob, enorm;
165 int i, j, iswap, jswap, nn, nuphill = 0;
166 gmx_rng_t rng;
167 t_mat *minimum;
168
169 if (m->n1 != m->nn)
170 {
171 fprintf(stderrstderr, "Can not do Monte Carlo optimization with a non-square matrix.\n");
172 return;
173 }
174 printf("\nDoing Monte Carlo optimization to find the smoothest trajectory\n");
175 printf("by reordering the frames to minimize the path between the two structures\n");
176 printf("that have the largest pairwise RMSD.\n");
177
178 iswap = jswap = -1;
179 enorm = m->mat[0][0];
180 for (i = 0; (i < m->n1); i++)
181 {
182 for (j = 0; (j < m->nn); j++)
183 {
184 if (m->mat[i][j] > enorm)
185 {
186 enorm = m->mat[i][j];
187 iswap = i;
188 jswap = j;
189 }
190 }
191 }
192 if ((iswap == -1) || (jswap == -1))
193 {
194 fprintf(stderrstderr, "Matrix contains identical values in all fields\n");
195 return;
196 }
197 swap_rows(m, 0, iswap);
198 swap_rows(m, m->n1-1, jswap);
199 emin = ecur = mat_energy(m);
200 printf("Largest distance %g between %d and %d. Energy: %g.\n",
201 enorm, iswap, jswap, emin);
202
203 rng = gmx_rng_init(seed);
204 nn = m->nn;
205
206 /* Initiate and store global minimum */
207 minimum = init_mat(nn, m->b1D);
208 minimum->nn = nn;
209 copy_t_mat(minimum, m);
210
211 if (NULL((void*)0) != conv)
212 {
213 fp = xvgropen(conv, "Convergence of the MC optimization",
214 "Energy", "Step", oenv);
215 }
216 for (i = 0; (i < maxiter); i++)
217 {
218 /* Generate new swapping candidates */
219 do
220 {
221 iswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
222 jswap = 1+(nn-2)*gmx_rng_uniform_real(rng);
223 }
224 while ((iswap == jswap) || (iswap >= nn-1) || (jswap >= nn-1));
225
226 /* Apply swap and compute energy */
227 swap_rows(m, iswap, jswap);
228 enext = mat_energy(m);
229
230 /* Compute probability */
231 prob = 0;
232 if ((enext < ecur) || (i < nrandom))
233 {
234 prob = 1;
235 if (enext < emin)
236 {
237 /* Store global minimum */
238 copy_t_mat(minimum, m);
239 emin = enext;
240 }
241 }
242 else if (kT > 0)
243 {
244 /* Try Monte Carlo step */
245 prob = exp(-(enext-ecur)/(enorm*kT));
246 }
247
248 if ((prob == 1) || (gmx_rng_uniform_real(rng) < prob))
249 {
250 if (enext > ecur)
251 {
252 nuphill++;
253 }
254
255 fprintf(log, "Iter: %d Swapped %4d and %4d (energy: %g prob: %g)\n",
256 i, iswap, jswap, enext, prob);
257 if (NULL((void*)0) != fp)
258 {
259 fprintf(fp, "%6d %10g\n", i, enext);
260 }
261 ecur = enext;
262 }
263 else
264 {
265 swap_rows(m, jswap, iswap);
266 }
267 }
268 fprintf(log, "%d uphill steps were taken during optimization\n",
269 nuphill);
270
271 /* Now swap the matrix to get it into global minimum mode */
272 copy_t_mat(m, minimum);
273
274 fprintf(log, "Global minimum energy %g\n", mat_energy(minimum));
275 fprintf(log, "Global minimum energy %g\n", mat_energy(m));
276 fprintf(log, "Swapped time and frame indices and RMSD to next neighbor:\n");
277 for (i = 0; (i < m->nn); i++)
278 {
279 fprintf(log, "%10g %5d %10g\n",
280 time[m->m_ind[i]],
281 m->m_ind[i],
282 (i < m->nn-1) ? m->mat[m->m_ind[i]][m->m_ind[i+1]] : 0);
283 }
284
285 if (NULL((void*)0) != fp)
286 {
287 fclose(fp);
288 }
289}
290
291static void calc_dist(int nind, rvec x[], real **d)
292{
293 int i, j;
294 real *xi;
295 rvec dx;
296
297 for (i = 0; (i < nind-1); i++)
298 {
299 xi = x[i];
300 for (j = i+1; (j < nind); j++)
301 {
302 /* Should use pbc_dx when analysing multiple molecueles,
303 * but the box is not stored for every frame.
304 */
305 rvec_sub(xi, x[j], dx);
306 d[i][j] = norm(dx);
307 }
308 }
309}
310
311static real rms_dist(int isize, real **d, real **d_r)
312{
313 int i, j;
314 real r, r2;
315
316 r2 = 0.0;
317 for (i = 0; (i < isize-1); i++)
318 {
319 for (j = i+1; (j < isize); j++)
320 {
321 r = d[i][j]-d_r[i][j];
322 r2 += r*r;
323 }
324 }
325 r2 /= (isize*(isize-1))/2;
326
327 return sqrt(r2);
328}
329
330static int rms_dist_comp(const void *a, const void *b)
331{
332 t_dist *da, *db;
333
334 da = (t_dist *)a;
335 db = (t_dist *)b;
336
337 if (da->dist - db->dist < 0)
338 {
339 return -1;
340 }
341 else if (da->dist - db->dist > 0)
342 {
343 return 1;
344 }
345 return 0;
346}
347
348static int clust_id_comp(const void *a, const void *b)
349{
350 t_clustid *da, *db;
351
352 da = (t_clustid *)a;
353 db = (t_clustid *)b;
354
355 return da->clust - db->clust;
356}
357
358static int nrnb_comp(const void *a, const void *b)
359{
360 t_nnb *da, *db;
361
362 da = (t_nnb *)a;
363 db = (t_nnb *)b;
364
365 /* return the b-a, we want highest first */
366 return db->nr - da->nr;
367}
368
369void gather(t_mat *m, real cutoff, t_clusters *clust)
370{
371 t_clustid *c;
372 t_dist *d;
373 int i, j, k, nn, cid, n1, diff;
374 gmx_bool bChange;
375
376 /* First we sort the entries in the RMSD matrix */
377 n1 = m->nn;
378 nn = ((n1-1)*n1)/2;
379 snew(d, nn)(d) = save_calloc("d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 379, (nn), sizeof(*(d)));
380 for (i = k = 0; (i < n1); i++)
381 {
382 for (j = i+1; (j < n1); j++, k++)
383 {
384 d[k].i = i;
385 d[k].j = j;
386 d[k].dist = m->mat[i][j];
387 }
388 }
389 if (k != nn)
390 {
391 gmx_incons("gather algortihm")_gmx_error("incons", "gather algortihm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 391);
392 }
393 qsort(d, nn, sizeof(d[0]), rms_dist_comp);
394
395 /* Now we make a cluster index for all of the conformations */
396 c = new_clustid(n1);
397
398 /* Now we check the closest structures, and equalize their cluster numbers */
399 fprintf(stderrstderr, "Linking structures ");
400 do
401 {
402 fprintf(stderrstderr, "*");
403 bChange = FALSE0;
404 for (k = 0; (k < nn) && (d[k].dist < cutoff); k++)
405 {
406 diff = c[d[k].j].clust - c[d[k].i].clust;
407 if (diff)
408 {
409 bChange = TRUE1;
410 if (diff > 0)
411 {
412 c[d[k].j].clust = c[d[k].i].clust;
413 }
414 else
415 {
416 c[d[k].i].clust = c[d[k].j].clust;
417 }
418 }
419 }
420 }
421 while (bChange);
422 fprintf(stderrstderr, "\nSorting and renumbering clusters\n");
423 /* Sort on cluster number */
424 qsort(c, n1, sizeof(c[0]), clust_id_comp);
425
426 /* Renumber clusters */
427 cid = 1;
428 for (k = 1; k < n1; k++)
429 {
430 if (c[k].clust != c[k-1].clust)
431 {
432 c[k-1].clust = cid;
433 cid++;
434 }
435 else
436 {
437 c[k-1].clust = cid;
438 }
439 }
440 c[k-1].clust = cid;
441 if (debug)
442 {
443 for (k = 0; (k < n1); k++)
444 {
445 fprintf(debug, "Cluster index for conformation %d: %d\n",
446 c[k].conf, c[k].clust);
447 }
448 }
449 clust->ncl = cid;
450 for (k = 0; k < n1; k++)
451 {
452 clust->cl[c[k].conf] = c[k].clust;
453 }
454
455 sfree(c)save_free("c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 455, (c));
456 sfree(d)save_free("d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 456, (d));
457}
458
459gmx_bool jp_same(int **nnb, int i, int j, int P)
460{
461 gmx_bool bIn;
462 int k, ii, jj, pp;
463
464 bIn = FALSE0;
465 for (k = 0; nnb[i][k] >= 0; k++)
466 {
467 bIn = bIn || (nnb[i][k] == j);
468 }
469 if (!bIn)
470 {
471 return FALSE0;
472 }
473
474 bIn = FALSE0;
475 for (k = 0; nnb[j][k] >= 0; k++)
476 {
477 bIn = bIn || (nnb[j][k] == i);
478 }
479 if (!bIn)
480 {
481 return FALSE0;
482 }
483
484 pp = 0;
485 for (ii = 0; nnb[i][ii] >= 0; ii++)
486 {
487 for (jj = 0; nnb[j][jj] >= 0; jj++)
488 {
489 if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
490 {
491 pp++;
492 }
493 }
494 }
495
496 return (pp >= P);
497}
498
499static void jarvis_patrick(int n1, real **mat, int M, int P,
500 real rmsdcut, t_clusters *clust)
501{
502 t_dist *row;
503 t_clustid *c;
504 int **nnb;
505 int i, j, k, cid, diff, max;
506 gmx_bool bChange;
507 real **mcpy = NULL((void*)0);
508
509 if (rmsdcut < 0)
510 {
511 rmsdcut = 10000;
512 }
513
514 /* First we sort the entries in the RMSD matrix row by row.
515 * This gives us the nearest neighbor list.
516 */
517 snew(nnb, n1)(nnb) = save_calloc("nnb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 517, (n1), sizeof(*(nnb)));
518 snew(row, n1)(row) = save_calloc("row", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 518, (n1), sizeof(*(row)));
519 for (i = 0; (i < n1); i++)
520 {
521 for (j = 0; (j < n1); j++)
522 {
523 row[j].j = j;
524 row[j].dist = mat[i][j];
525 }
526 qsort(row, n1, sizeof(row[0]), rms_dist_comp);
527 if (M > 0)
528 {
529 /* Put the M nearest neighbors in the list */
530 snew(nnb[i], M+1)(nnb[i]) = save_calloc("nnb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 530, (M+1), sizeof(*(nnb[i])));
531 for (j = k = 0; (k < M) && (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
532 {
533 if (row[j].j != i)
534 {
535 nnb[i][k] = row[j].j;
536 k++;
537 }
538 }
539 nnb[i][k] = -1;
540 }
541 else
542 {
543 /* Put all neighbors nearer than rmsdcut in the list */
544 max = 0;
545 k = 0;
546 for (j = 0; (j < n1) && (mat[i][row[j].j] < rmsdcut); j++)
547 {
548 if (row[j].j != i)
549 {
550 if (k >= max)
551 {
552 max += 10;
553 srenew(nnb[i], max)(nnb[i]) = save_realloc("nnb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 553, (nnb[i]), (max), sizeof(*(nnb[i])));
554 }
555 nnb[i][k] = row[j].j;
556 k++;
557 }
558 }
559 if (k == max)
560 {
561 srenew(nnb[i], max+1)(nnb[i]) = save_realloc("nnb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 561, (nnb[i]), (max+1), sizeof(*(nnb[i])));
562 }
563 nnb[i][k] = -1;
564 }
565 }
566 sfree(row)save_free("row", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 566, (row));
567 if (debug)
568 {
569 fprintf(debug, "Nearest neighborlist. M = %d, P = %d\n", M, P);
570 for (i = 0; (i < n1); i++)
571 {
572 fprintf(debug, "i:%5d nbs:", i);
573 for (j = 0; nnb[i][j] >= 0; j++)
574 {
575 fprintf(debug, "%5d[%5.3f]", nnb[i][j], mat[i][nnb[i][j]]);
576 }
577 fprintf(debug, "\n");
578 }
579 }
580
581 c = new_clustid(n1);
582 fprintf(stderrstderr, "Linking structures ");
583 /* Use mcpy for temporary storage of booleans */
584 mcpy = mk_matrix(n1, n1, FALSE0);
585 for (i = 0; i < n1; i++)
586 {
587 for (j = i+1; j < n1; j++)
588 {
589 mcpy[i][j] = jp_same(nnb, i, j, P);
590 }
591 }
592 do
593 {
594 fprintf(stderrstderr, "*");
595 bChange = FALSE0;
596 for (i = 0; i < n1; i++)
597 {
598 for (j = i+1; j < n1; j++)
599 {
600 if (mcpy[i][j])
601 {
602 diff = c[j].clust - c[i].clust;
603 if (diff)
604 {
605 bChange = TRUE1;
606 if (diff > 0)
607 {
608 c[j].clust = c[i].clust;
609 }
610 else
611 {
612 c[i].clust = c[j].clust;
613 }
614 }
615 }
616 }
617 }
618 }
619 while (bChange);
620
621 fprintf(stderrstderr, "\nSorting and renumbering clusters\n");
622 /* Sort on cluster number */
623 qsort(c, n1, sizeof(c[0]), clust_id_comp);
624
625 /* Renumber clusters */
626 cid = 1;
627 for (k = 1; k < n1; k++)
628 {
629 if (c[k].clust != c[k-1].clust)
630 {
631 c[k-1].clust = cid;
632 cid++;
633 }
634 else
635 {
636 c[k-1].clust = cid;
637 }
638 }
639 c[k-1].clust = cid;
640 clust->ncl = cid;
641 for (k = 0; k < n1; k++)
642 {
643 clust->cl[c[k].conf] = c[k].clust;
644 }
645 if (debug)
646 {
647 for (k = 0; (k < n1); k++)
648 {
649 fprintf(debug, "Cluster index for conformation %d: %d\n",
650 c[k].conf, c[k].clust);
651 }
652 }
653
654/* Again, I don't see the point in this... (AF) */
655/* for(i=0; (i<n1); i++) { */
656/* for(j=0; (j<n1); j++) */
657/* mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
658/* } */
659/* for(i=0; (i<n1); i++) { */
660/* for(j=0; (j<n1); j++) */
661/* mat[i][j] = mcpy[i][j]; */
662/* } */
663 done_matrix(n1, &mcpy);
664
665 sfree(c)save_free("c", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 665, (c));
666 for (i = 0; (i < n1); i++)
667 {
668 sfree(nnb[i])save_free("nnb[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 668, (nnb[i]));
669 }
670 sfree(nnb)save_free("nnb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 670, (nnb));
671}
672
673static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
674{
675 int i, j;
676
677 /* dump neighbor list */
678 fprintf(fp, "%s", title);
679 for (i = 0; (i < n1); i++)
680 {
681 fprintf(fp, "i:%5d #:%5d nbs:", i, nnb[i].nr);
682 for (j = 0; j < nnb[i].nr; j++)
683 {
684 fprintf(fp, "%5d", nnb[i].nb[j]);
685 }
686 fprintf(fp, "\n");
687 }
688}
689
690static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
691{
692 t_dist *row;
693 t_nnb *nnb;
694 int i, j, k, j1, max;
695
696 /* Put all neighbors nearer than rmsdcut in the list */
697 fprintf(stderrstderr, "Making list of neighbors within cutoff ");
698 snew(nnb, n1)(nnb) = save_calloc("nnb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 698, (n1), sizeof(*(nnb)));
699 snew(row, n1)(row) = save_calloc("row", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 699, (n1), sizeof(*(row)));
700 for (i = 0; (i < n1); i++)
701 {
702 max = 0;
703 k = 0;
704 /* put all neighbors within cut-off in list */
705 for (j = 0; j < n1; j++)
706 {
707 if (mat[i][j] < rmsdcut)
708 {
709 if (k >= max)
710 {
711 max += 10;
712 srenew(nnb[i].nb, max)(nnb[i].nb) = save_realloc("nnb[i].nb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 712, (nnb[i].nb), (max), sizeof(*(nnb[i].nb)));
713 }
714 nnb[i].nb[k] = j;
715 k++;
716 }
717 }
718 /* store nr of neighbors, we'll need that */
719 nnb[i].nr = k;
720 if (i%(1+n1/100) == 0)
721 {
722 fprintf(stderrstderr, "%3d%%\b\b\b\b", (i*100+1)/n1);
723 }
724 }
725 fprintf(stderrstderr, "%3d%%\n", 100);
726 sfree(row)save_free("row", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 726, (row));
727
728 /* sort neighbor list on number of neighbors, largest first */
729 qsort(nnb, n1, sizeof(nnb[0]), nrnb_comp);
730
731 if (debug)
732 {
733 dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
734 }
735
736 /* turn first structure with all its neighbors (largest) into cluster
737 remove them from pool of structures and repeat for all remaining */
738 fprintf(stderrstderr, "Finding clusters %4d", 0);
739 /* cluster id's start at 1: */
740 k = 1;
741 while (nnb[0].nr)
742 {
743 /* set cluster id (k) for first item in neighborlist */
744 for (j = 0; j < nnb[0].nr; j++)
745 {
746 clust->cl[nnb[0].nb[j]] = k;
747 }
748 /* mark as done */
749 nnb[0].nr = 0;
750 sfree(nnb[0].nb)save_free("nnb[0].nb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 750, (nnb[0].nb));
751
752 /* adjust number of neighbors for others, taking removals into account: */
753 for (i = 1; i < n1 && nnb[i].nr; i++)
754 {
755 j1 = 0;
756 for (j = 0; j < nnb[i].nr; j++)
757 {
758 /* if this neighbor wasn't removed */
759 if (clust->cl[nnb[i].nb[j]] == 0)
760 {
761 /* shift the rest (j1<=j) */
762 nnb[i].nb[j1] = nnb[i].nb[j];
763 /* next */
764 j1++;
765 }
766 }
767 /* now j1 is the new number of neighbors */
768 nnb[i].nr = j1;
769 }
770 /* sort again on nnb[].nr, because we have new # neighbors: */
771 /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
772 qsort(nnb, i, sizeof(nnb[0]), nrnb_comp);
773
774 fprintf(stderrstderr, "\b\b\b\b%4d", k);
775 /* new cluster id */
776 k++;
777 }
778 fprintf(stderrstderr, "\n");
779 sfree(nnb)save_free("nnb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 779, (nnb));
780 if (debug)
781 {
782 fprintf(debug, "Clusters (%d):\n", k);
783 for (i = 0; i < n1; i++)
784 {
785 fprintf(debug, " %3d", clust->cl[i]);
786 }
787 fprintf(debug, "\n");
788 }
789
790 clust->ncl = k-1;
791}
792
793rvec **read_whole_trj(const char *fn, int isize, atom_id index[], int skip,
794 int *nframe, real **time, const output_env_t oenv, gmx_bool bPBC, gmx_rmpbc_t gpbc)
795{
796 rvec **xx, *x;
797 matrix box;
798 real t;
799 int i, i0, j, max_nf;
800 int natom;
801 t_trxstatus *status;
802
803
804 max_nf = 0;
805 xx = NULL((void*)0);
806 *time = NULL((void*)0);
807 natom = read_first_x(oenv, &status, fn, &t, &x, box);
808 i = 0;
809 i0 = 0;
810 do
811 {
812 if (bPBC)
813 {
814 gmx_rmpbc(gpbc, natom, box, x);
815 }
816 if (i0 >= max_nf)
817 {
818 max_nf += 10;
819 srenew(xx, max_nf)(xx) = save_realloc("xx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 819, (xx), (max_nf), sizeof(*(xx)));
820 srenew(*time, max_nf)(*time) = save_realloc("*time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 820, (*time), (max_nf), sizeof(*(*time)));
821 }
822 if ((i % skip) == 0)
823 {
824 snew(xx[i0], isize)(xx[i0]) = save_calloc("xx[i0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 824, (isize), sizeof(*(xx[i0])));
825 /* Store only the interesting atoms */
826 for (j = 0; (j < isize); j++)
827 {
828 copy_rvec(x[index[j]], xx[i0][j]);
829 }
830 (*time)[i0] = t;
831 i0++;
832 }
833 i++;
834 }
835 while (read_next_x(oenv, status, &t, x, box));
836 fprintf(stderrstderr, "Allocated %lu bytes for frames\n",
837 (unsigned long) (max_nf*isize*sizeof(**xx)));
838 fprintf(stderrstderr, "Read %d frames from trajectory %s\n", i0, fn);
839 *nframe = i0;
840 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 840, (x));
841
842 return xx;
843}
844
845static int plot_clusters(int nf, real **mat, t_clusters *clust,
846 int minstruct)
847{
848 int i, j, ncluster, ci;
849 int *cl_id, *nstruct, *strind;
850
851 snew(cl_id, nf)(cl_id) = save_calloc("cl_id", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 851, (nf), sizeof(*(cl_id)));
852 snew(nstruct, nf)(nstruct) = save_calloc("nstruct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 852, (nf), sizeof(*(nstruct)));
853 snew(strind, nf)(strind) = save_calloc("strind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 853, (nf), sizeof(*(strind)));
854 for (i = 0; i < nf; i++)
855 {
856 strind[i] = 0;
857 cl_id[i] = clust->cl[i];
858 nstruct[cl_id[i]]++;
859 }
860 ncluster = 0;
861 for (i = 0; i < nf; i++)
862 {
863 if (nstruct[i] >= minstruct)
864 {
865 ncluster++;
866 for (j = 0; (j < nf); j++)
867 {
868 if (cl_id[j] == i)
869 {
870 strind[j] = ncluster;
871 }
872 }
873 }
874 }
875 ncluster++;
876 fprintf(stderrstderr, "There are %d clusters with at least %d conformations\n",
877 ncluster, minstruct);
878
879 for (i = 0; (i < nf); i++)
880 {
881 ci = cl_id[i];
882 for (j = 0; j < i; j++)
883 {
884 if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct))
885 {
886 /* color different clusters with different colors, as long as
887 we don't run out of colors */
888 mat[i][j] = strind[i];
889 }
890 else
891 {
892 mat[i][j] = 0;
893 }
894 }
895 }
896 sfree(strind)save_free("strind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 896, (strind));
897 sfree(nstruct)save_free("nstruct", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 897, (nstruct));
898 sfree(cl_id)save_free("cl_id", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 898, (cl_id));
899
900 return ncluster;
901}
902
903static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
904{
905 int i, j, v;
906
907 for (i = 0; i < nf; i++)
908 {
909 for (j = 0; j < i; j++)
910 {
911 if (clust->cl[i] == clust->cl[j])
912 {
913 mat[i][j] = val;
914 }
915 else
916 {
917 mat[i][j] = 0;
918 }
919 }
920 }
921}
922
923static char *parse_filename(const char *fn, int maxnr)
924{
925 int i;
926 char *fnout;
927 const char *ext;
928 char buf[STRLEN4096];
929
930 if (strchr(fn, '%')(__extension__ (__builtin_constant_p ('%') && !__builtin_constant_p
(fn) && ('%') == '\0' ? (char *) __rawmemchr (fn, '%'
) : __builtin_strchr (fn, '%'))))
931 {
932 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 932, "will not number filename %s containing '%c'", fn, '%');
933 }
934 /* number of digits needed in numbering */
935 i = (int)(log(maxnr)/log(10)) + 1;
936 /* split fn and ext */
937 ext = strrchr(fn, '.');
938 if (!ext)
939 {
940 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 940, "cannot separate extension in filename %s", fn);
941 }
942 ext++;
943 /* insert e.g. '%03d' between fn and ext */
944 sprintf(buf, "%s%%0%dd.%s", fn, i, ext);
945 snew(fnout, strlen(buf)+1)(fnout) = save_calloc("fnout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 945, (strlen(buf)+1), sizeof(*(fnout)));
946 strcpy(fnout, buf);
947
948 return fnout;
949}
950
951static void ana_trans(t_clusters *clust, int nf,
952 const char *transfn, const char *ntransfn, FILE *log,
953 t_rgb rlo, t_rgb rhi, const output_env_t oenv)
954{
955 FILE *fp;
956 real **trans, *axis;
957 int *ntrans;
958 int i, ntranst, maxtrans;
959 char buf[STRLEN4096];
960
961 snew(ntrans, clust->ncl)(ntrans) = save_calloc("ntrans", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 961, (clust->ncl), sizeof(*(ntrans)));
962 snew(trans, clust->ncl)(trans) = save_calloc("trans", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 962, (clust->ncl), sizeof(*(trans)));
963 snew(axis, clust->ncl)(axis) = save_calloc("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 963, (clust->ncl), sizeof(*(axis)));
964 for (i = 0; i < clust->ncl; i++)
965 {
966 axis[i] = i+1;
967 snew(trans[i], clust->ncl)(trans[i]) = save_calloc("trans[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 967, (clust->ncl), sizeof(*(trans[i])));
968 }
969 ntranst = 0;
970 maxtrans = 0;
971 for (i = 1; i < nf; i++)
972 {
973 if (clust->cl[i] != clust->cl[i-1])
974 {
975 ntranst++;
976 ntrans[clust->cl[i-1]-1]++;
977 ntrans[clust->cl[i]-1]++;
978 trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
979 maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1])(((maxtrans) > (trans[clust->cl[i]-1][clust->cl[i-1]
-1])) ? (maxtrans) : (trans[clust->cl[i]-1][clust->cl[i
-1]-1]) );
980 }
981 }
982 ffprintf_dd(stderrstderr, log, buf, "Counted %d transitions in total, "
983 "max %d between two specific clusters\n", ntranst, maxtrans);
984 if (transfn)
985 {
986 fp = gmx_ffopen(transfn, "w");
987 i = min(maxtrans+1, 80)(((maxtrans+1) < (80)) ? (maxtrans+1) : (80) );
988 write_xpm(fp, 0, "Cluster Transitions", "# transitions",
989 "from cluster", "to cluster",
990 clust->ncl, clust->ncl, axis, axis, trans,
991 0, maxtrans, rlo, rhi, &i);
992 gmx_ffclose(fp);
993 }
994 if (ntransfn)
995 {
996 fp = xvgropen(ntransfn, "Cluster Transitions", "Cluster #", "# transitions",
997 oenv);
998 for (i = 0; i < clust->ncl; i++)
999 {
1000 fprintf(fp, "%5d %5d\n", i+1, ntrans[i]);
1001 }
1002 gmx_ffclose(fp);
1003 }
1004 sfree(ntrans)save_free("ntrans", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1004, (ntrans));
1005 for (i = 0; i < clust->ncl; i++)
1006 {
1007 sfree(trans[i])save_free("trans[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1007, (trans[i]));
1008 }
1009 sfree(trans)save_free("trans", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1009, (trans));
1010 sfree(axis)save_free("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1010, (axis));
1011}
1012
1013static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
1014 int natom, t_atoms *atoms, rvec *xtps,
1015 real *mass, rvec **xx, real *time,
1016 int ifsize, atom_id *fitidx,
1017 int iosize, atom_id *outidx,
1018 const char *trxfn, const char *sizefn,
1019 const char *transfn, const char *ntransfn,
1020 const char *clustidfn, gmx_bool bAverage,
1021 int write_ncl, int write_nst, real rmsmin,
1022 gmx_bool bFit, FILE *log, t_rgb rlo, t_rgb rhi,
1023 const output_env_t oenv)
1024{
1025 FILE *fp = NULL((void*)0);
1026 char buf[STRLEN4096], buf1[40], buf2[40], buf3[40], *trxsfn;
1027 t_trxstatus *trxout = NULL((void*)0);
1028 t_trxstatus *trxsout = NULL((void*)0);
1029 int i, i1, cl, nstr, *structure, first = 0, midstr;
1030 gmx_bool *bWrite = NULL((void*)0);
1031 real r, clrmsd, midrmsd;
1032 rvec *xav = NULL((void*)0);
1033 matrix zerobox;
1034
1035 clear_mat(zerobox);
1036
1037 ffprintf_d(stderrstderr, log, buf, "\nFound %d clusters\n\n", clust->ncl);
1038 trxsfn = NULL((void*)0);
1039 if (trxfn)
1040 {
1041 /* do we write all structures? */
1042 if (write_ncl)
1043 {
1044 trxsfn = parse_filename(trxfn, max(write_ncl, clust->ncl)(((write_ncl) > (clust->ncl)) ? (write_ncl) : (clust->
ncl) ));
1045 snew(bWrite, nf)(bWrite) = save_calloc("bWrite", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1045, (nf), sizeof(*(bWrite)));
1046 }
1047 ffprintf_ss(stderrstderr, log, buf, "Writing %s structure for each cluster to %s\n",
1048 bAverage ? "average" : "middle", trxfn);
1049 if (write_ncl)
1050 {
1051 /* find out what we want to tell the user:
1052 Writing [all structures|structures with rmsd > %g] for
1053 {all|first %d} clusters {with more than %d structures} to %s */
1054 if (rmsmin > 0.0)
1055 {
1056 sprintf(buf1, "structures with rmsd > %g", rmsmin);
1057 }
1058 else
1059 {
1060 sprintf(buf1, "all structures");
1061 }
1062 buf2[0] = buf3[0] = '\0';
1063 if (write_ncl >= clust->ncl)
1064 {
1065 if (write_nst == 0)
1066 {
1067 sprintf(buf2, "all ");
1068 }
1069 }
1070 else
1071 {
1072 sprintf(buf2, "the first %d ", write_ncl);
1073 }
1074 if (write_nst)
1075 {
1076 sprintf(buf3, " with more than %d structures", write_nst);
1077 }
1078 sprintf(buf, "Writing %s for %sclusters%s to %s\n", buf1, buf2, buf3, trxsfn);
1079 ffprintf(stderrstderr, log, buf);
1080 }
1081
1082 /* Prepare a reference structure for the orientation of the clusters */
1083 if (bFit)
1084 {
1085 reset_x(ifsize, fitidx, natom, NULL((void*)0), xtps, mass);
1086 }
1087 trxout = open_trx(trxfn, "w");
1088 /* Calculate the average structure in each cluster, *
1089 * all structures are fitted to the first struture of the cluster */
1090 snew(xav, natom)(xav) = save_calloc("xav", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1090, (natom), sizeof(*(xav)));
1091 }
1092
1093 if (transfn || ntransfn)
1094 {
1095 ana_trans(clust, nf, transfn, ntransfn, log, rlo, rhi, oenv);
1096 }
1097
1098 if (clustidfn)
1099 {
1100 fp = xvgropen(clustidfn, "Clusters", output_env_get_xvgr_tlabel(oenv), "Cluster #", oenv);
1101 fprintf(fp, "@ s0 symbol 2\n");
1102 fprintf(fp, "@ s0 symbol size 0.2\n");
1103 fprintf(fp, "@ s0 linestyle 0\n");
1104 for (i = 0; i < nf; i++)
1105 {
1106 fprintf(fp, "%8g %8d\n", time[i], clust->cl[i]);
1107 }
1108 gmx_ffclose(fp);
1109 }
1110 if (sizefn)
1111 {
1112 fp = xvgropen(sizefn, "Cluster Sizes", "Cluster #", "# Structures", oenv);
1113 fprintf(fp, "@g%d type %s\n", 0, "bar");
1114 }
1115 snew(structure, nf)(structure) = save_calloc("structure", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1115, (nf), sizeof(*(structure)));
1116 fprintf(log, "\n%3s | %3s %4s | %6s %4s | cluster members\n",
1117 "cl.", "#st", "rmsd", "middle", "rmsd");
1118 for (cl = 1; cl <= clust->ncl; cl++)
1119 {
1120 /* prepare structures (fit, middle, average) */
1121 if (xav)
1122 {
1123 for (i = 0; i < natom; i++)
1124 {
1125 clear_rvec(xav[i]);
1126 }
1127 }
1128 nstr = 0;
1129 for (i1 = 0; i1 < nf; i1++)
1130 {
1131 if (clust->cl[i1] == cl)
1132 {
1133 structure[nstr] = i1;
1134 nstr++;
1135 if (trxfn && (bAverage || write_ncl) )
1136 {
1137 if (bFit)
1138 {
1139 reset_x(ifsize, fitidx, natom, NULL((void*)0), xx[i1], mass);
1140 }
1141 if (nstr == 1)
1142 {
1143 first = i1;
1144 }
1145 else if (bFit)
1146 {
1147 do_fit(natom, mass, xx[first], xx[i1]);
1148 }
1149 if (xav)
1150 {
1151 for (i = 0; i < natom; i++)
1152 {
1153 rvec_inc(xav[i], xx[i1][i]);
1154 }
1155 }
1156 }
1157 }
1158 }
1159 if (sizefn)
1160 {
1161 fprintf(fp, "%8d %8d\n", cl, nstr);
1162 }
1163 clrmsd = 0;
1164 midstr = 0;
1165 midrmsd = 10000;
1166 for (i1 = 0; i1 < nstr; i1++)
1167 {
1168 r = 0;
1169 if (nstr > 1)
1170 {
1171 for (i = 0; i < nstr; i++)
1172 {
1173 if (i < i1)
1174 {
1175 r += rmsd[structure[i]][structure[i1]];
1176 }
1177 else
1178 {
1179 r += rmsd[structure[i1]][structure[i]];
1180 }
1181 }
1182 r /= (nstr - 1);
1183 }
1184 if (r < midrmsd)
1185 {
1186 midstr = structure[i1];
1187 midrmsd = r;
1188 }
1189 clrmsd += r;
1190 }
1191 clrmsd /= nstr;
1192
1193 /* dump cluster info to logfile */
1194 if (nstr > 1)
1195 {
1196 sprintf(buf1, "%6.3f", clrmsd);
1197 if (buf1[0] == '0')
1198 {
1199 buf1[0] = ' ';
1200 }
1201 sprintf(buf2, "%5.3f", midrmsd);
1202 if (buf2[0] == '0')
1203 {
1204 buf2[0] = ' ';
1205 }
1206 }
1207 else
1208 {
1209 sprintf(buf1, "%5s", "");
1210 sprintf(buf2, "%5s", "");
1211 }
1212 fprintf(log, "%3d | %3d %s | %6g%s |", cl, nstr, buf1, time[midstr], buf2);
1213 for (i = 0; i < nstr; i++)
1214 {
1215 if ((i % 7 == 0) && i)
1216 {
1217 sprintf(buf, "\n%3s | %3s %4s | %6s %4s |", "", "", "", "", "");
1218 }
1219 else
1220 {
1221 buf[0] = '\0';
1222 }
1223 i1 = structure[i];
1224 fprintf(log, "%s %6g", buf, time[i1]);
1225 }
1226 fprintf(log, "\n");
1227
1228 /* write structures to trajectory file(s) */
1229 if (trxfn)
1230 {
1231 if (write_ncl)
1232 {
1233 for (i = 0; i < nstr; i++)
1234 {
1235 bWrite[i] = FALSE0;
1236 }
1237 }
1238 if (cl < write_ncl+1 && nstr > write_nst)
1239 {
1240 /* Dump all structures for this cluster */
1241 /* generate numbered filename (there is a %d in trxfn!) */
1242 sprintf(buf, trxsfn, cl);
1243 trxsout = open_trx(buf, "w");
1244 for (i = 0; i < nstr; i++)
1245 {
1246 bWrite[i] = TRUE1;
1247 if (rmsmin > 0.0)
1248 {
1249 for (i1 = 0; i1 < i && bWrite[i]; i1++)
1250 {
1251 if (bWrite[i1])
1252 {
1253 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
1254 }
1255 }
1256 }
1257 if (bWrite[i])
1258 {
1259 write_trx(trxsout, iosize, outidx, atoms, i, time[structure[i]], zerobox,
1260 xx[structure[i]], NULL((void*)0), NULL((void*)0));
1261 }
1262 }
1263 close_trx(trxsout);
1264 }
1265 /* Dump the average structure for this cluster */
1266 if (bAverage)
1267 {
1268 for (i = 0; i < natom; i++)
1269 {
1270 svmul(1.0/nstr, xav[i], xav[i]);
1271 }
1272 }
1273 else
1274 {
1275 for (i = 0; i < natom; i++)
1276 {
1277 copy_rvec(xx[midstr][i], xav[i]);
1278 }
1279 if (bFit)
1280 {
1281 reset_x(ifsize, fitidx, natom, NULL((void*)0), xav, mass);
1282 }
1283 }
1284 if (bFit)
1285 {
1286 do_fit(natom, mass, xtps, xav);
1287 }
1288 r = cl;
Value stored to 'r' is never read
1289 write_trx(trxout, iosize, outidx, atoms, cl, time[midstr], zerobox, xav, NULL((void*)0), NULL((void*)0));
1290 }
1291 }
1292 /* clean up */
1293 if (trxfn)
1294 {
1295 close_trx(trxout);
1296 sfree(xav)save_free("xav", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1296, (xav));
1297 if (write_ncl)
1298 {
1299 sfree(bWrite)save_free("bWrite", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1299, (bWrite));
1300 }
1301 }
1302 sfree(structure)save_free("structure", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1302, (structure));
1303 if (trxsfn)
1304 {
1305 sfree(trxsfn)save_free("trxsfn", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1305, (trxsfn));
1306 }
1307}
1308
1309static void convert_mat(t_matrix *mat, t_mat *rms)
1310{
1311 int i, j;
1312
1313 rms->n1 = mat->nx;
1314 matrix2real(mat, rms->mat);
1315 /* free input xpm matrix data */
1316 for (i = 0; i < mat->nx; i++)
1317 {
1318 sfree(mat->matrix[i])save_free("mat->matrix[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1318, (mat->matrix[i]));
1319 }
1320 sfree(mat->matrix)save_free("mat->matrix", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1320, (mat->matrix));
1321
1322 for (i = 0; i < mat->nx; i++)
1323 {
1324 for (j = i; j < mat->nx; j++)
1325 {
1326 rms->sumrms += rms->mat[i][j];
1327 rms->maxrms = max(rms->maxrms, rms->mat[i][j])(((rms->maxrms) > (rms->mat[i][j])) ? (rms->maxrms
) : (rms->mat[i][j]) );
1328 if (j != i)
1329 {
1330 rms->minrms = min(rms->minrms, rms->mat[i][j])(((rms->minrms) < (rms->mat[i][j])) ? (rms->minrms
) : (rms->mat[i][j]) );
1331 }
1332 }
1333 }
1334 rms->nn = mat->nx;
1335}
1336
1337int gmx_cluster(int argc, char *argv[])
1338{
1339 const char *desc[] = {
1340 "[THISMODULE] can cluster structures using several different methods.",
1341 "Distances between structures can be determined from a trajectory",
1342 "or read from an [TT].xpm[tt] matrix file with the [TT]-dm[tt] option.",
1343 "RMS deviation after fitting or RMS deviation of atom-pair distances",
1344 "can be used to define the distance between structures.[PAR]",
1345
1346 "single linkage: add a structure to a cluster when its distance to any",
1347 "element of the cluster is less than [TT]cutoff[tt].[PAR]",
1348
1349 "Jarvis Patrick: add a structure to a cluster when this structure",
1350 "and a structure in the cluster have each other as neighbors and",
1351 "they have a least [TT]P[tt] neighbors in common. The neighbors",
1352 "of a structure are the M closest structures or all structures within",
1353 "[TT]cutoff[tt].[PAR]",
1354
1355 "Monte Carlo: reorder the RMSD matrix using Monte Carlo such that",
1356 "the order of the frames is using the smallest possible increments.",
1357 "With this it is possible to make a smooth animation going from one",
1358 "structure to another with the largest possible (e.g.) RMSD between",
1359 "them, however the intermediate steps should be as small as possible.",
1360 "Applications could be to visualize a potential of mean force",
1361 "ensemble of simulations or a pulling simulation. Obviously the user",
1362 "has to prepare the trajectory well (e.g. by not superimposing frames).",
1363 "The final result can be inspect visually by looking at the matrix",
1364 "[TT].xpm[tt] file, which should vary smoothly from bottom to top.[PAR]",
1365
1366 "diagonalization: diagonalize the RMSD matrix.[PAR]",
1367
1368 "gromos: use algorithm as described in Daura [IT]et al.[it]",
1369 "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
1370 "Count number of neighbors using cut-off, take structure with",
1371 "largest number of neighbors with all its neighbors as cluster",
1372 "and eliminate it from the pool of clusters. Repeat for remaining",
1373 "structures in pool.[PAR]",
1374
1375 "When the clustering algorithm assigns each structure to exactly one",
1376 "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
1377 "file is supplied, the structure with",
1378 "the smallest average distance to the others or the average structure",
1379 "or all structures for each cluster will be written to a trajectory",
1380 "file. When writing all structures, separate numbered files are made",
1381 "for each cluster.[PAR]",
1382
1383 "Two output files are always written:[BR]",
1384 "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
1385 "and a graphical depiction of the clusters in the lower right half",
1386 "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
1387 "when two structures are in the same cluster.",
1388 "When [TT]-minstruct[tt] > 1 different colors will be used for each",
1389 "cluster.[BR]",
1390 "[TT]-g[tt] writes information on the options used and a detailed list",
1391 "of all clusters and their members.[PAR]",
1392
1393 "Additionally, a number of optional output files can be written:[BR]",
1394 "[TT]-dist[tt] writes the RMSD distribution.[BR]",
1395 "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
1396 "diagonalization.[BR]",
1397 "[TT]-sz[tt] writes the cluster sizes.[BR]",
1398 "[TT]-tr[tt] writes a matrix of the number transitions between",
1399 "cluster pairs.[BR]",
1400 "[TT]-ntr[tt] writes the total number of transitions to or from",
1401 "each cluster.[BR]",
1402 "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1403 "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1404 "structure of each cluster or writes numbered files with cluster members",
1405 "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1406 "[TT]-nst[tt] and [TT]-rmsmin[tt]). The center of a cluster is the",
1407 "structure with the smallest average RMSD from all other structures",
1408 "of the cluster.[BR]",
1409 };
1410
1411 FILE *fp, *log;
1412 int nf, i, i1, i2, j;
1413 gmx_int64_t nrms = 0;
1414
1415 matrix box;
1416 rvec *xtps, *usextps, *x1, **xx = NULL((void*)0);
1417 const char *fn, *trx_out_fn;
1418 t_clusters clust;
1419 t_mat *rms, *orig = NULL((void*)0);
1420 real *eigenvalues;
1421 t_topology top;
1422 int ePBC;
1423 t_atoms useatoms;
1424 t_matrix *readmat = NULL((void*)0);
1425 real *eigenvectors;
1426
1427 int isize = 0, ifsize = 0, iosize = 0;
1428 atom_id *index = NULL((void*)0), *fitidx, *outidx;
1429 char *grpname;
1430 real rmsd, **d1, **d2, *time = NULL((void*)0), time_invfac, *mass = NULL((void*)0);
1431 char buf[STRLEN4096], buf1[80], title[STRLEN4096];
1432 gmx_bool bAnalyze, bUseRmsdCut, bJP_RMSD = FALSE0, bReadMat, bReadTraj, bPBC = TRUE1;
1433
1434 int method, ncluster = 0;
1435 static const char *methodname[] = {
1436 NULL((void*)0), "linkage", "jarvis-patrick", "monte-carlo",
1437 "diagonalization", "gromos", NULL((void*)0)
1438 };
1439 enum {
1440 m_null, m_linkage, m_jarvis_patrick,
1441 m_monte_carlo, m_diagonalize, m_gromos, m_nr
1442 };
1443 /* Set colors for plotting: white = zero RMS, black = maximum */
1444 static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1445 static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1446 static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1447 static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1448 static int nlevels = 40, skip = 1;
1449 static real scalemax = -1.0, rmsdcut = 0.1, rmsmin = 0.0;
1450 gmx_bool bRMSdist = FALSE0, bBinary = FALSE0, bAverage = FALSE0, bFit = TRUE1;
1451 static int niter = 10000, nrandom = 0, seed = 1993, write_ncl = 0, write_nst = 1, minstruct = 1;
1452 static real kT = 1e-3;
1453 static int M = 10, P = 3;
1454 output_env_t oenv;
1455 gmx_rmpbc_t gpbc = NULL((void*)0);
1456
1457 t_pargs pa[] = {
1458 { "-dista", FALSE0, etBOOL, {&bRMSdist},
1459 "Use RMSD of distances instead of RMS deviation" },
1460 { "-nlevels", FALSE0, etINT, {&nlevels},
1461 "Discretize RMSD matrix in this number of levels" },
1462 { "-cutoff", FALSE0, etREAL, {&rmsdcut},
1463 "RMSD cut-off (nm) for two structures to be neighbor" },
1464 { "-fit", FALSE0, etBOOL, {&bFit},
1465 "Use least squares fitting before RMSD calculation" },
1466 { "-max", FALSE0, etREAL, {&scalemax},
1467 "Maximum level in RMSD matrix" },
1468 { "-skip", FALSE0, etINT, {&skip},
1469 "Only analyze every nr-th frame" },
1470 { "-av", FALSE0, etBOOL, {&bAverage},
1471 "Write average iso middle structure for each cluster" },
1472 { "-wcl", FALSE0, etINT, {&write_ncl},
1473 "Write the structures for this number of clusters to numbered files" },
1474 { "-nst", FALSE0, etINT, {&write_nst},
1475 "Only write all structures if more than this number of structures per cluster" },
1476 { "-rmsmin", FALSE0, etREAL, {&rmsmin},
1477 "minimum rms difference with rest of cluster for writing structures" },
1478 { "-method", FALSE0, etENUM, {methodname},
1479 "Method for cluster determination" },
1480 { "-minstruct", FALSE0, etINT, {&minstruct},
1481 "Minimum number of structures in cluster for coloring in the [TT].xpm[tt] file" },
1482 { "-binary", FALSE0, etBOOL, {&bBinary},
1483 "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1484 "is given by [TT]-cutoff[tt]" },
1485 { "-M", FALSE0, etINT, {&M},
1486 "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1487 "0 is use cutoff" },
1488 { "-P", FALSE0, etINT, {&P},
1489 "Number of identical nearest neighbors required to form a cluster" },
1490 { "-seed", FALSE0, etINT, {&seed},
1491 "Random number seed for Monte Carlo clustering algorithm: <= 0 means generate" },
1492 { "-niter", FALSE0, etINT, {&niter},
1493 "Number of iterations for MC" },
1494 { "-nrandom", FALSE0, etINT, {&nrandom},
1495 "The first iterations for MC may be done complete random, to shuffle the frames" },
1496 { "-kT", FALSE0, etREAL, {&kT},
1497 "Boltzmann weighting factor for Monte Carlo optimization "
1498 "(zero turns off uphill steps)" },
1499 { "-pbc", FALSE0, etBOOL,
1500 { &bPBC }, "PBC check" }
1501 };
1502 t_filenm fnm[] = {
1503 { efTRX, "-f", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
1504 { efTPS, "-s", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
1505 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
1506 { efXPM, "-dm", "rmsd", ffOPTRD(1<<1 | 1<<3) },
1507 { efXPM, "-om", "rmsd-raw", ffWRITE1<<2 },
1508 { efXPM, "-o", "rmsd-clust", ffWRITE1<<2 },
1509 { efLOG, "-g", "cluster", ffWRITE1<<2 },
1510 { efXVG, "-dist", "rmsd-dist", ffOPTWR(1<<2| 1<<3) },
1511 { efXVG, "-ev", "rmsd-eig", ffOPTWR(1<<2| 1<<3) },
1512 { efXVG, "-conv", "mc-conv", ffOPTWR(1<<2| 1<<3) },
1513 { efXVG, "-sz", "clust-size", ffOPTWR(1<<2| 1<<3)},
1514 { efXPM, "-tr", "clust-trans", ffOPTWR(1<<2| 1<<3)},
1515 { efXVG, "-ntr", "clust-trans", ffOPTWR(1<<2| 1<<3)},
1516 { efXVG, "-clid", "clust-id.xvg", ffOPTWR(1<<2| 1<<3)},
1517 { efTRX, "-cl", "clusters.pdb", ffOPTWR(1<<2| 1<<3) }
1518 };
1519#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
1520
1521 if (!parse_common_args(&argc, argv,
1522 PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13),
1523 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0),
1524 &oenv))
1525 {
1526 return 0;
1527 }
1528
1529 /* parse options */
1530 bReadMat = opt2bSet("-dm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1531 bReadTraj = opt2bSet("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || !bReadMat;
1532 if (opt2parg_bSet("-av", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ||
1533 opt2parg_bSet("-wcl", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ||
1534 opt2parg_bSet("-nst", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ||
1535 opt2parg_bSet("-rmsmin", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) ||
1536 opt2bSet("-cl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) )
1537 {
1538 trx_out_fn = opt2fn("-cl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1539 }
1540 else
1541 {
1542 trx_out_fn = NULL((void*)0);
1543 }
1544 if (bReadMat && output_env_get_time_factor(oenv) != 1)
1545 {
1546 fprintf(stderrstderr,
1547 "\nWarning: assuming the time unit in %s is %s\n",
1548 opt2fn("-dm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), output_env_get_time_unit(oenv));
1549 }
1550 if (trx_out_fn && !bReadTraj)
1551 {
1552 fprintf(stderrstderr, "\nWarning: "
1553 "cannot write cluster structures without reading trajectory\n"
1554 " ignoring option -cl %s\n", trx_out_fn);
1555 }
1556
1557 method = 1;
1558 while (method < m_nr && gmx_strcasecmp(methodname[0], methodname[method]) != 0)
1559 {
1560 method++;
1561 }
1562 if (method == m_nr)
1563 {
1564 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1564, "Invalid method");
1565 }
1566
1567 bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1568 method == m_gromos );
1569
1570 /* Open log file */
1571 log = ftp2FILE(efLOG, NFILE, fnm, "w")gmx_ffopen(ftp2fn(efLOG, ((int)(sizeof(fnm)/sizeof((fnm)[0]))
), fnm), "w");
1572
1573 fprintf(stderrstderr, "Using %s method for clustering\n", methodname[0]);
1574 fprintf(log, "Using %s method for clustering\n", methodname[0]);
1575
1576 /* check input and write parameters to log file */
1577 bUseRmsdCut = FALSE0;
1578 if (method == m_jarvis_patrick)
1579 {
1580 bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa);
1581 if ((M < 0) || (M == 1))
1582 {
1583 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1583, "M (%d) must be 0 or larger than 1", M);
1584 }
1585 if (M < 2)
1586 {
1587 sprintf(buf1, "Will use P=%d and RMSD cutoff (%g)", P, rmsdcut);
1588 bUseRmsdCut = TRUE1;
1589 }
1590 else
1591 {
1592 if (P >= M)
1593 {
1594 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1594, "Number of neighbors required (P) must be less than M");
1595 }
1596 if (bJP_RMSD)
1597 {
1598 sprintf(buf1, "Will use P=%d, M=%d and RMSD cutoff (%g)", P, M, rmsdcut);
1599 bUseRmsdCut = TRUE1;
1600 }
1601 else
1602 {
1603 sprintf(buf1, "Will use P=%d, M=%d", P, M);
1604 }
1605 }
1606 ffprintf_s(stderrstderr, log, buf, "%s for determining the neighbors\n\n", buf1);
1607 }
1608 else /* method != m_jarvis */
1609 {
1610 bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1611 }
1612 if (bUseRmsdCut && method != m_jarvis_patrick)
1613 {
1614 fprintf(log, "Using RMSD cutoff %g nm\n", rmsdcut);
1615 }
1616 if (method == m_monte_carlo)
1617 {
1618 fprintf(log, "Using %d iterations\n", niter);
1619 }
1620
1621 if (skip < 1)
1622 {
1623 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1623, "skip (%d) should be >= 1", skip);
1624 }
1625
1626 /* get input */
1627 if (bReadTraj)
1628 {
1629 /* don't read mass-database as masses (and top) are not used */
1630 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, &top, &ePBC, &xtps, NULL((void*)0), box,
1631 TRUE1);
1632 if (bPBC)
1633 {
1634 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
1635 }
1636
1637 fprintf(stderrstderr, "\nSelect group for least squares fit%s:\n",
1638 bReadMat ? "" : " and RMSD calculation");
1639 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1640 1, &ifsize, &fitidx, &grpname);
1641 if (trx_out_fn)
1642 {
1643 fprintf(stderrstderr, "\nSelect group for output:\n");
1644 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1645 1, &iosize, &outidx, &grpname);
1646 /* merge and convert both index groups: */
1647 /* first copy outidx to index. let outidx refer to elements in index */
1648 snew(index, iosize)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1648, (iosize), sizeof(*(index)));
1649 isize = iosize;
1650 for (i = 0; i < iosize; i++)
1651 {
1652 index[i] = outidx[i];
1653 outidx[i] = i;
1654 }
1655 /* now lookup elements from fitidx in index, add them if necessary
1656 and also let fitidx refer to elements in index */
1657 for (i = 0; i < ifsize; i++)
1658 {
1659 j = 0;
1660 while (j < isize && index[j] != fitidx[i])
1661 {
1662 j++;
1663 }
1664 if (j >= isize)
1665 {
1666 /* slow this way, but doesn't matter much */
1667 isize++;
1668 srenew(index, isize)(index) = save_realloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1668, (index), (isize), sizeof(*(index)));
1669 }
1670 index[j] = fitidx[i];
1671 fitidx[i] = j;
1672 }
1673 }
1674 else /* !trx_out_fn */
1675 {
1676 isize = ifsize;
1677 snew(index, isize)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1677, (isize), sizeof(*(index)));
1678 for (i = 0; i < ifsize; i++)
1679 {
1680 index[i] = fitidx[i];
1681 fitidx[i] = i;
1682 }
1683 }
1684 }
1685
1686 if (bReadTraj)
1687 {
1688 /* Loop over first coordinate file */
1689 fn = opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
1690
1691 xx = read_whole_trj(fn, isize, index, skip, &nf, &time, oenv, bPBC, gpbc);
1692 output_env_conv_times(oenv, nf, time);
1693 if (!bRMSdist || bAnalyze)
1694 {
1695 /* Center all frames on zero */
1696 snew(mass, isize)(mass) = save_calloc("mass", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1696, (isize), sizeof(*(mass)));
1697 for (i = 0; i < ifsize; i++)
1698 {
1699 mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1700 }
1701 if (bFit)
1702 {
1703 for (i = 0; i < nf; i++)
1704 {
1705 reset_x(ifsize, fitidx, isize, NULL((void*)0), xx[i], mass);
1706 }
1707 }
1708 }
1709 if (bPBC)
1710 {
1711 gmx_rmpbc_done(gpbc);
1712 }
1713 }
1714
1715 if (bReadMat)
1716 {
1717 fprintf(stderrstderr, "Reading rms distance matrix ");
1718 read_xpm_matrix(opt2fn("-dm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &readmat);
1719 fprintf(stderrstderr, "\n");
1720 if (readmat[0].nx != readmat[0].ny)
1721 {
1722 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1722, "Matrix (%dx%d) is not square",
1723 readmat[0].nx, readmat[0].ny);
1724 }
1725 if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1726 {
1727 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1727, "Matrix size (%dx%d) does not match the number of "
1728 "frames (%d)", readmat[0].nx, readmat[0].ny, nf);
1729 }
1730
1731 nf = readmat[0].nx;
1732 sfree(time)save_free("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1732, (time));
1733 time = readmat[0].axis_x;
1734 time_invfac = output_env_get_time_invfactor(oenv);
1735 for (i = 0; i < nf; i++)
1736 {
1737 time[i] *= time_invfac;
1738 }
1739
1740 rms = init_mat(readmat[0].nx, method == m_diagonalize);
1741 convert_mat(&(readmat[0]), rms);
1742
1743 nlevels = readmat[0].nmap;
1744 }
1745 else /* !bReadMat */
1746 {
1747 rms = init_mat(nf, method == m_diagonalize);
1748 nrms = ((gmx_int64_t)nf*((gmx_int64_t)nf-1))/2;
1749 if (!bRMSdist)
1750 {
1751 fprintf(stderrstderr, "Computing %dx%d RMS deviation matrix\n", nf, nf);
1752 /* Initialize work array */
1753 snew(x1, isize)(x1) = save_calloc("x1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1753, (isize), sizeof(*(x1)));
1754 for (i1 = 0; i1 < nf; i1++)
1755 {
1756 for (i2 = i1+1; i2 < nf; i2++)
1757 {
1758 for (i = 0; i < isize; i++)
1759 {
1760 copy_rvec(xx[i1][i], x1[i]);
1761 }
1762 if (bFit)
1763 {
1764 do_fit(isize, mass, xx[i2], x1);
1765 }
1766 rmsd = rmsdev(isize, mass, xx[i2], x1);
1767 set_mat_entry(rms, i1, i2, rmsd);
1768 }
1769 nrms -= (gmx_int64_t) (nf-i1-1);
1770 fprintf(stderrstderr, "\r# RMSD calculations left: " "%"GMX_PRId64"l" "d" " ", nrms);
1771 }
1772 sfree(x1)save_free("x1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1772, (x1));
1773 }
1774 else /* bRMSdist */
1775 {
1776 fprintf(stderrstderr, "Computing %dx%d RMS distance deviation matrix\n", nf, nf);
1777
1778 /* Initiate work arrays */
1779 snew(d1, isize)(d1) = save_calloc("d1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1779, (isize), sizeof(*(d1)));
1780 snew(d2, isize)(d2) = save_calloc("d2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1780, (isize), sizeof(*(d2)));
1781 for (i = 0; (i < isize); i++)
1782 {
1783 snew(d1[i], isize)(d1[i]) = save_calloc("d1[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1783, (isize), sizeof(*(d1[i])));
1784 snew(d2[i], isize)(d2[i]) = save_calloc("d2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1784, (isize), sizeof(*(d2[i])));
1785 }
1786 for (i1 = 0; i1 < nf; i1++)
1787 {
1788 calc_dist(isize, xx[i1], d1);
1789 for (i2 = i1+1; (i2 < nf); i2++)
1790 {
1791 calc_dist(isize, xx[i2], d2);
1792 set_mat_entry(rms, i1, i2, rms_dist(isize, d1, d2));
1793 }
1794 nrms -= (nf-i1-1);
1795 fprintf(stderrstderr, "\r# RMSD calculations left: " "%"GMX_PRId64"l" "d" " ", nrms);
1796 }
1797 /* Clean up work arrays */
1798 for (i = 0; (i < isize); i++)
1799 {
1800 sfree(d1[i])save_free("d1[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1800, (d1[i]));
1801 sfree(d2[i])save_free("d2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1801, (d2[i]));
1802 }
1803 sfree(d1)save_free("d1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1803, (d1));
1804 sfree(d2)save_free("d2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1804, (d2));
1805 }
1806 fprintf(stderrstderr, "\n\n");
1807 }
1808 ffprintf_gg(stderrstderr, log, buf, "The RMSD ranges from %g to %g nm\n",
1809 rms->minrms, rms->maxrms);
1810 ffprintf_g(stderrstderr, log, buf, "Average RMSD is %g\n", 2*rms->sumrms/(nf*(nf-1)));
1811 ffprintf_d(stderrstderr, log, buf, "Number of structures for matrix %d\n", nf);
1812 ffprintf_g(stderrstderr, log, buf, "Energy of the matrix is %g.\n", mat_energy(rms));
1813 if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1814 {
1815 fprintf(stderrstderr, "WARNING: rmsd cutoff %g is outside range of rmsd values "
1816 "%g to %g\n", rmsdcut, rms->minrms, rms->maxrms);
1817 }
1818 if (bAnalyze && (rmsmin < rms->minrms) )
1819 {
1820 fprintf(stderrstderr, "WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1821 rmsmin, rms->minrms);
1822 }
1823 if (bAnalyze && (rmsmin > rmsdcut) )
1824 {
1825 fprintf(stderrstderr, "WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1826 rmsmin, rmsdcut);
1827 }
1828
1829 /* Plot the rmsd distribution */
1830 rmsd_distribution(opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), rms, oenv);
1831
1832 if (bBinary)
1833 {
1834 for (i1 = 0; (i1 < nf); i1++)
1835 {
1836 for (i2 = 0; (i2 < nf); i2++)
1837 {
1838 if (rms->mat[i1][i2] < rmsdcut)
1839 {
1840 rms->mat[i1][i2] = 0;
1841 }
1842 else
1843 {
1844 rms->mat[i1][i2] = 1;
1845 }
1846 }
1847 }
1848 }
1849
1850 snew(clust.cl, nf)(clust.cl) = save_calloc("clust.cl", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1850, (nf), sizeof(*(clust.cl)));
1851 switch (method)
1852 {
1853 case m_linkage:
1854 /* Now sort the matrix and write it out again */
1855 gather(rms, rmsdcut, &clust);
1856 break;
1857 case m_diagonalize:
1858 /* Do a diagonalization */
1859 snew(eigenvalues, nf)(eigenvalues) = save_calloc("eigenvalues", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1859, (nf), sizeof(*(eigenvalues)));
1860 snew(eigenvectors, nf*nf)(eigenvectors) = save_calloc("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1860, (nf*nf), sizeof(*(eigenvectors)));
1861 memcpy(eigenvectors, rms->mat[0], nf*nf*sizeof(real));
1862 eigensolver(eigenvectors, nf, 0, nf, eigenvalues, rms->mat[0]);
1863 sfree(eigenvectors)save_free("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1863, (eigenvectors));
1864
1865 fp = xvgropen(opt2fn("-ev", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "RMSD matrix Eigenvalues",
1866 "Eigenvector index", "Eigenvalues (nm\\S2\\N)", oenv);
1867 for (i = 0; (i < nf); i++)
1868 {
1869 fprintf(fp, "%10d %10g\n", i, eigenvalues[i]);
1870 }
1871 gmx_ffclose(fp);
1872 break;
1873 case m_monte_carlo:
1874 orig = init_mat(rms->nn, FALSE0);
1875 orig->nn = rms->nn;
1876 copy_t_mat(orig, rms);
1877 mc_optimize(log, rms, time, niter, nrandom, seed, kT,
1878 opt2fn_null("-conv", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), oenv);
1879 break;
1880 case m_jarvis_patrick:
1881 jarvis_patrick(rms->nn, rms->mat, M, P, bJP_RMSD ? rmsdcut : -1, &clust);
1882 break;
1883 case m_gromos:
1884 gromos(rms->nn, rms->mat, rmsdcut, &clust);
1885 break;
1886 default:
1887 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1887, "DEATH HORROR unknown method \"%s\"", methodname[0]);
1888 }
1889
1890 if (method == m_monte_carlo || method == m_diagonalize)
1891 {
1892 fprintf(stderrstderr, "Energy of the matrix after clustering is %g.\n",
1893 mat_energy(rms));
1894 }
1895
1896 if (bAnalyze)
1897 {
1898 if (minstruct > 1)
1899 {
1900 ncluster = plot_clusters(nf, rms->mat, &clust, minstruct);
1901 }
1902 else
1903 {
1904 mark_clusters(nf, rms->mat, rms->maxrms, &clust);
1905 }
1906 init_t_atoms(&useatoms, isize, FALSE0);
1907 snew(usextps, isize)(usextps) = save_calloc("usextps", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1907, (isize), sizeof(*(usextps)));
1908 useatoms.resinfo = top.atoms.resinfo;
1909 for (i = 0; i < isize; i++)
1910 {
1911 useatoms.atomname[i] = top.atoms.atomname[index[i]];
1912 useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1913 useatoms.nres = max(useatoms.nres, useatoms.atom[i].resind+1)(((useatoms.nres) > (useatoms.atom[i].resind+1)) ? (useatoms
.nres) : (useatoms.atom[i].resind+1) );
1914 copy_rvec(xtps[index[i]], usextps[i]);
1915 }
1916 useatoms.nr = isize;
1917 analyze_clusters(nf, &clust, rms->mat, isize, &useatoms, usextps, mass, xx, time,
1918 ifsize, fitidx, iosize, outidx,
1919 bReadTraj ? trx_out_fn : NULL((void*)0),
1920 opt2fn_null("-sz", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1921 opt2fn_null("-tr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1922 opt2fn_null("-ntr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1923 opt2fn_null("-clid", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
1924 bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1925 rlo_bot, rhi_bot, oenv);
1926 }
1927 gmx_ffclose(log);
1928
1929 if (bBinary && !bAnalyze)
1930 {
1931 /* Make the clustering visible */
1932 for (i2 = 0; (i2 < nf); i2++)
1933 {
1934 for (i1 = i2+1; (i1 < nf); i1++)
1935 {
1936 if (rms->mat[i1][i2])
1937 {
1938 rms->mat[i1][i2] = rms->maxrms;
1939 }
1940 }
1941 }
1942 }
1943
1944 fp = opt2FILE("-o", NFILE, fnm, "w")gmx_ffopen(opt2fn("-o", ((int)(sizeof(fnm)/sizeof((fnm)[0])))
, fnm), "w");
1945 fprintf(stderrstderr, "Writing rms distance/clustering matrix ");
1946 if (bReadMat)
1947 {
1948 write_xpm(fp, 0, readmat[0].title, readmat[0].legend, readmat[0].label_x,
1949 readmat[0].label_y, nf, nf, readmat[0].axis_x, readmat[0].axis_y,
1950 rms->mat, 0.0, rms->maxrms, rlo_top, rhi_top, &nlevels);
1951 }
1952 else
1953 {
1954 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1955 sprintf(title, "RMS%sDeviation / Cluster Index",
1956 bRMSdist ? " Distance " : " ");
1957 if (minstruct > 1)
1958 {
1959 write_xpm_split(fp, 0, title, "RMSD (nm)", buf, buf,
1960 nf, nf, time, time, rms->mat, 0.0, rms->maxrms, &nlevels,
1961 rlo_top, rhi_top, 0.0, (real) ncluster,
1962 &ncluster, TRUE1, rlo_bot, rhi_bot);
1963 }
1964 else
1965 {
1966 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1967 nf, nf, time, time, rms->mat, 0.0, rms->maxrms,
1968 rlo_top, rhi_top, &nlevels);
1969 }
1970 }
1971 fprintf(stderrstderr, "\n");
1972 gmx_ffclose(fp);
1973 if (NULL((void*)0) != orig)
1974 {
1975 fp = opt2FILE("-om", NFILE, fnm, "w")gmx_ffopen(opt2fn("-om", ((int)(sizeof(fnm)/sizeof((fnm)[0]))
), fnm), "w");
1976 sprintf(buf, "Time (%s)", output_env_get_time_unit(oenv));
1977 sprintf(title, "RMS%sDeviation", bRMSdist ? " Distance " : " ");
1978 write_xpm(fp, 0, title, "RMSD (nm)", buf, buf,
1979 nf, nf, time, time, orig->mat, 0.0, orig->maxrms,
1980 rlo_top, rhi_top, &nlevels);
1981 gmx_ffclose(fp);
1982 done_mat(&orig);
1983 sfree(orig)save_free("orig", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_cluster.c"
, 1983, (orig));
1984 }
1985 /* now show what we've done */
1986 do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1987 do_view(oenv, opt2fn_null("-sz", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1988 if (method == m_diagonalize)
1989 {
1990 do_view(oenv, opt2fn_null("-ev", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1991 }
1992 do_view(oenv, opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1993 if (bAnalyze)
1994 {
1995 do_view(oenv, opt2fn_null("-tr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1996 do_view(oenv, opt2fn_null("-ntr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1997 do_view(oenv, opt2fn_null("-clid", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
1998 }
1999 do_view(oenv, opt2fn_null("-conv", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
2000
2001 return 0;
2002}