File: | gromacs/gmxana/gmx_cluster.c |
Location: | line 1908, column 26 |
Description: | Assigned value is garbage or undefined |
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) */ | |||
71 | static gmx_inlineinline | |||
72 | void 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 */ | |||
79 | static gmx_inlineinline | |||
80 | void 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 */ | |||
86 | static gmx_inlineinline | |||
87 | void 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 */ | |||
94 | static gmx_inlineinline | |||
95 | void 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 */ | |||
102 | static gmx_inlineinline | |||
103 | void 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 */ | |||
110 | static gmx_inlineinline | |||
111 | void 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 */ | |||
118 | static gmx_inlineinline | |||
119 | void 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 */ | |||
126 | static gmx_inlineinline | |||
127 | void 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 | ||||
133 | typedef struct { | |||
134 | int ncl; | |||
135 | int *cl; | |||
136 | } t_clusters; | |||
137 | ||||
138 | typedef struct { | |||
139 | int nr; | |||
140 | int *nb; | |||
141 | } t_nnb; | |||
142 | ||||
143 | void pr_energy(FILE *fp, real e) | |||
144 | { | |||
145 | fprintf(fp, "Energy: %8.4f\n", e); | |||
146 | } | |||
147 | ||||
148 | void 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 | ||||
158 | void 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 | ||||
291 | static 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 | ||||
311 | static 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 | ||||
330 | static 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 | ||||
348 | static 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 | ||||
358 | static 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 | ||||
369 | void 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 | ||||
459 | gmx_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 | ||||
499 | static 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 | ||||
673 | static 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 | ||||
690 | static 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 | ||||
793 | rvec **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 | ||||
845 | static 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 | ||||
903 | static 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 | ||||
923 | static 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 | ||||
951 | static 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 | ||||
1013 | static 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; | |||
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 | ||||
1309 | static 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 | ||||
1337 | int 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 | } |