File: | gromacs/gmxana/gmx_anaeig.c |
Location: | line 95, column 13 |
Description: | Value stored to 'w' is never read |
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 "gromacs/commandline/pargs.h" |
46 | #include "typedefs.h" |
47 | #include "gromacs/utility/smalloc.h" |
48 | #include "macros.h" |
49 | #include "gromacs/utility/fatalerror.h" |
50 | #include "gromacs/math/vec.h" |
51 | #include "pbc.h" |
52 | #include "gromacs/utility/futil.h" |
53 | #include "index.h" |
54 | #include "gromacs/fileio/pdbio.h" |
55 | #include "gromacs/fileio/confio.h" |
56 | #include "gromacs/fileio/tpxio.h" |
57 | #include "gromacs/fileio/trxio.h" |
58 | #include "gromacs/fileio/matio.h" |
59 | #include "mshift.h" |
60 | #include "gromacs/fileio/xvgr.h" |
61 | #include "viewit.h" |
62 | #include "rmpbc.h" |
63 | #include "txtdump.h" |
64 | #include "eigio.h" |
65 | #include "physics.h" |
66 | #include "gmx_ana.h" |
67 | |
68 | #include "gromacs/math/do_fit.h" |
69 | |
70 | static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip) |
71 | { |
72 | int i; |
73 | double hwkT, w, dS, S = 0; |
74 | double hbar, lambda; |
75 | |
76 | hbar = PLANCK1(6.6262e-34)/(2*M_PI3.14159265358979323846); |
77 | for (i = 0; (i < n-nskip); i++) |
78 | { |
79 | if (eigval[i] > 0) |
80 | { |
81 | lambda = eigval[i]*AMU(1.6605402e-27); |
82 | w = sqrt(BOLTZMANN(1.380658e-23)*temp/lambda)/NANO(1e-9); |
83 | hwkT = (hbar*w)/(BOLTZMANN(1.380658e-23)*temp); |
84 | dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT))); |
85 | S += dS; |
86 | if (debug) |
87 | { |
88 | fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n", |
89 | i, w, lambda, hwkT, dS); |
90 | } |
91 | } |
92 | else |
93 | { |
94 | fprintf(stderrstderr, "eigval[%d] = %g\n", i, eigval[i]); |
95 | w = 0; |
Value stored to 'w' is never read | |
96 | } |
97 | } |
98 | fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n", |
99 | S*RGAS((1.380658e-23)*(6.0221367e23))); |
100 | } |
101 | |
102 | static void calc_entropy_schlitter(FILE *fp, int n, int nskip, |
103 | real *eigval, real temp) |
104 | { |
105 | double dd, deter; |
106 | int *indx; |
107 | int i, j, k, m; |
108 | char buf[256]; |
109 | double hbar, kt, kteh, S; |
110 | |
111 | hbar = PLANCK1(6.6262e-34)/(2*M_PI3.14159265358979323846); |
112 | kt = BOLTZMANN(1.380658e-23)*temp; |
113 | kteh = kt*exp(2.0)/(hbar*hbar)*AMU(1.6605402e-27)*(NANO(1e-9)*NANO(1e-9)); |
114 | if (debug) |
115 | { |
116 | fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh); |
117 | } |
118 | |
119 | deter = 0; |
120 | for (i = 0; (i < n-nskip); i++) |
121 | { |
122 | dd = 1+kteh*eigval[i]; |
123 | deter += log(dd); |
124 | } |
125 | S = 0.5*RGAS((1.380658e-23)*(6.0221367e23))*deter; |
126 | |
127 | fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S); |
128 | } |
129 | |
130 | const char *proj_unit; |
131 | |
132 | static real tick_spacing(real range, int minticks) |
133 | { |
134 | real sp; |
135 | |
136 | if (range <= 0) |
137 | { |
138 | return 1.0; |
139 | } |
140 | |
141 | sp = 0.2*exp(log(10)*ceil(log(range)/log(10))); |
142 | while (range/sp < minticks-1) |
143 | { |
144 | sp = sp/2; |
145 | } |
146 | |
147 | return sp; |
148 | } |
149 | |
150 | static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph, |
151 | const char *title, const char *subtitle, |
152 | const char *xlabel, const char **ylabel, |
153 | int n, real *x, real **y, real ***sy, |
154 | real scale_x, gmx_bool bZero, gmx_bool bSplit, |
155 | const output_env_t oenv) |
156 | { |
157 | FILE *out; |
158 | int g, s, i; |
159 | real min, max, xsp, ysp; |
160 | |
161 | out = gmx_ffopen(file, "w"); |
162 | if (output_env_get_xvg_format(oenv) == exvgXMGRACE) |
163 | { |
164 | fprintf(out, "@ autoscale onread none\n"); |
165 | } |
166 | for (g = 0; g < ngraphs; g++) |
167 | { |
168 | if (y) |
169 | { |
170 | min = y[g][0]; |
171 | max = y[g][0]; |
172 | for (i = 0; i < n; i++) |
173 | { |
174 | if (y[g][i] < min) |
175 | { |
176 | min = y[g][i]; |
177 | } |
178 | if (y[g][i] > max) |
179 | { |
180 | max = y[g][i]; |
181 | } |
182 | } |
183 | } |
184 | else |
185 | { |
186 | min = sy[g][0][0]; |
187 | max = sy[g][0][0]; |
188 | for (s = 0; s < nsetspergraph; s++) |
189 | { |
190 | for (i = 0; i < n; i++) |
191 | { |
192 | if (sy[g][s][i] < min) |
193 | { |
194 | min = sy[g][s][i]; |
195 | } |
196 | if (sy[g][s][i] > max) |
197 | { |
198 | max = sy[g][s][i]; |
199 | } |
200 | } |
201 | } |
202 | } |
203 | if (bZero) |
204 | { |
205 | min = 0; |
206 | } |
207 | else |
208 | { |
209 | min = min-0.1*(max-min); |
210 | } |
211 | max = max+0.1*(max-min); |
212 | xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4); |
213 | ysp = tick_spacing(max-min, 3); |
214 | if (output_env_get_print_xvgr_codes(oenv)) |
215 | { |
216 | fprintf(out, "@ with g%d\n@ g%d on\n", g, g); |
217 | if (g == 0) |
218 | { |
219 | fprintf(out, "@ title \"%s\"\n", title); |
220 | if (subtitle) |
221 | { |
222 | fprintf(out, "@ subtitle \"%s\"\n", subtitle); |
223 | } |
224 | } |
225 | if (g == ngraphs-1) |
226 | { |
227 | fprintf(out, "@ xaxis label \"%s\"\n", xlabel); |
228 | } |
229 | else |
230 | { |
231 | fprintf(out, "@ xaxis ticklabel off\n"); |
232 | } |
233 | if (n > 1) |
234 | { |
235 | fprintf(out, "@ world xmin %g\n", x[0]*scale_x); |
236 | fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x); |
237 | fprintf(out, "@ world ymin %g\n", min); |
238 | fprintf(out, "@ world ymax %g\n", max); |
239 | } |
240 | fprintf(out, "@ view xmin 0.15\n"); |
241 | fprintf(out, "@ view xmax 0.85\n"); |
242 | fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs); |
243 | fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs); |
244 | fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]); |
245 | fprintf(out, "@ xaxis tick major %g\n", xsp); |
246 | fprintf(out, "@ xaxis tick minor %g\n", xsp/2); |
247 | fprintf(out, "@ xaxis ticklabel start type spec\n"); |
248 | fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp); |
249 | fprintf(out, "@ yaxis tick major %g\n", ysp); |
250 | fprintf(out, "@ yaxis tick minor %g\n", ysp/2); |
251 | fprintf(out, "@ yaxis ticklabel start type spec\n"); |
252 | fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp); |
253 | if ((min < 0) && (max > 0)) |
254 | { |
255 | fprintf(out, "@ zeroxaxis bar on\n"); |
256 | fprintf(out, "@ zeroxaxis bar linestyle 3\n"); |
257 | } |
258 | } |
259 | for (s = 0; s < nsetspergraph; s++) |
260 | { |
261 | for (i = 0; i < n; i++) |
262 | { |
263 | if (bSplit && i > 0 && abs(x[i]) < 1e-5) |
264 | { |
265 | if (output_env_get_print_xvgr_codes(oenv)) |
266 | { |
267 | fprintf(out, "&\n"); |
268 | } |
269 | } |
270 | fprintf(out, "%10.4f %10.5f\n", |
271 | x[i]*scale_x, y ? y[g][i] : sy[g][s][i]); |
272 | } |
273 | if (output_env_get_print_xvgr_codes(oenv)) |
274 | { |
275 | fprintf(out, "&\n"); |
276 | } |
277 | } |
278 | } |
279 | gmx_ffclose(out); |
280 | } |
281 | |
282 | static void |
283 | compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2, |
284 | real *eigval1, int neig1, real *eigval2, int neig2) |
285 | { |
286 | int n, nrow; |
287 | int i, j, k; |
288 | double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip; |
289 | |
290 | n = min(n1, n2)(((n1) < (n2)) ? (n1) : (n2) ); |
291 | |
292 | n = min(n, min(neig1, neig2))(((n) < ((((neig1) < (neig2)) ? (neig1) : (neig2) ))) ? (n) : ((((neig1) < (neig2)) ? (neig1) : (neig2) )) ); |
293 | fprintf(stdoutstdout, "Will compare the covariance matrices using %d dimensions\n", n); |
294 | |
295 | sum1 = 0; |
296 | for (i = 0; i < n; i++) |
297 | { |
298 | if (eigval1[i] < 0) |
299 | { |
300 | eigval1[i] = 0; |
301 | } |
302 | sum1 += eigval1[i]; |
303 | eigval1[i] = sqrt(eigval1[i]); |
304 | } |
305 | trace1 = sum1; |
306 | for (i = n; i < neig1; i++) |
307 | { |
308 | trace1 += eigval1[i]; |
309 | } |
310 | sum2 = 0; |
311 | for (i = 0; i < n; i++) |
312 | { |
313 | if (eigval2[i] < 0) |
314 | { |
315 | eigval2[i] = 0; |
316 | } |
317 | sum2 += eigval2[i]; |
318 | eigval2[i] = sqrt(eigval2[i]); |
319 | } |
320 | trace2 = sum2; |
321 | for (i = n; i < neig2; i++) |
322 | { |
323 | trace2 += eigval2[i]; |
324 | } |
325 | |
326 | fprintf(stdoutstdout, "Trace of the two matrices: %g and %g\n", sum1, sum2); |
327 | if (neig1 != n || neig2 != n) |
328 | { |
329 | fprintf(stdoutstdout, "this is %d%% and %d%% of the total trace\n", |
330 | (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5)); |
331 | } |
332 | fprintf(stdoutstdout, "Square root of the traces: %g and %g\n", |
333 | sqrt(sum1), sqrt(sum2)); |
334 | |
335 | sab = 0; |
336 | for (i = 0; i < n; i++) |
337 | { |
338 | tmp = 0; |
339 | for (j = 0; j < n; j++) |
340 | { |
341 | ip = 0; |
342 | for (k = 0; k < natoms; k++) |
343 | { |
344 | ip += iprod(eigvec1[i][k], eigvec2[j][k]); |
345 | } |
346 | tmp += eigval2[j]*ip*ip; |
347 | } |
348 | sab += eigval1[i]*tmp; |
349 | } |
350 | |
351 | samsb2 = sum1+sum2-2*sab; |
352 | if (samsb2 < 0) |
353 | { |
354 | samsb2 = 0; |
355 | } |
356 | |
357 | fprintf(stdoutstdout, "The overlap of the covariance matrices:\n"); |
358 | fprintf(stdoutstdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2))); |
359 | tmp = 1-sab/sqrt(sum1*sum2); |
360 | if (tmp < 0) |
361 | { |
362 | tmp = 0; |
363 | } |
364 | fprintf(stdoutstdout, " shape: %.3f\n", 1-sqrt(tmp)); |
365 | } |
366 | |
367 | |
368 | static void inprod_matrix(const char *matfile, int natoms, |
369 | int nvec1, int *eignr1, rvec **eigvec1, |
370 | int nvec2, int *eignr2, rvec **eigvec2, |
371 | gmx_bool bSelect, int noutvec, int *outvec) |
372 | { |
373 | FILE *out; |
374 | real **mat; |
375 | int i, x1, y1, x, y, nlevels; |
376 | int nx, ny; |
377 | real inp, *t_x, *t_y, max; |
378 | t_rgb rlo, rhi; |
379 | |
380 | snew(t_y, nvec2)(t_y) = save_calloc("t_y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 380, (nvec2), sizeof(*(t_y))); |
381 | if (bSelect) |
382 | { |
383 | nx = noutvec; |
384 | ny = 0; |
385 | for (y1 = 0; y1 < nx; y1++) |
386 | { |
387 | if (outvec[y1] < nvec2) |
388 | { |
389 | t_y[ny] = eignr2[outvec[y1]]+1; |
390 | ny++; |
391 | } |
392 | } |
393 | } |
394 | else |
395 | { |
396 | nx = nvec1; |
397 | ny = nvec2; |
398 | for (y = 0; y < ny; y++) |
399 | { |
400 | t_y[y] = eignr2[y]+1; |
401 | } |
402 | } |
403 | |
404 | fprintf(stderrstderr, "Calculating inner-product matrix of %dx%d eigenvectors\n", |
405 | nx, nvec2); |
406 | |
407 | snew(mat, nx)(mat) = save_calloc("mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 407, (nx), sizeof(*(mat))); |
408 | snew(t_x, nx)(t_x) = save_calloc("t_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 408, (nx), sizeof(*(t_x))); |
409 | max = 0; |
410 | for (x1 = 0; x1 < nx; x1++) |
411 | { |
412 | snew(mat[x1], ny)(mat[x1]) = save_calloc("mat[x1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 412, (ny), sizeof(*(mat[x1]))); |
413 | if (bSelect) |
414 | { |
415 | x = outvec[x1]; |
416 | } |
417 | else |
418 | { |
419 | x = x1; |
420 | } |
421 | t_x[x1] = eignr1[x]+1; |
422 | fprintf(stderrstderr, " %d", eignr1[x]+1); |
423 | for (y1 = 0; y1 < ny; y1++) |
424 | { |
425 | inp = 0; |
426 | if (bSelect) |
427 | { |
428 | while (outvec[y1] >= nvec2) |
429 | { |
430 | y1++; |
431 | } |
432 | y = outvec[y1]; |
433 | } |
434 | else |
435 | { |
436 | y = y1; |
437 | } |
438 | for (i = 0; i < natoms; i++) |
439 | { |
440 | inp += iprod(eigvec1[x][i], eigvec2[y][i]); |
441 | } |
442 | mat[x1][y1] = fabs(inp); |
443 | if (mat[x1][y1] > max) |
444 | { |
445 | max = mat[x1][y1]; |
446 | } |
447 | } |
448 | } |
449 | fprintf(stderrstderr, "\n"); |
450 | rlo.r = 1; rlo.g = 1; rlo.b = 1; |
451 | rhi.r = 0; rhi.g = 0; rhi.b = 0; |
452 | nlevels = 41; |
453 | out = gmx_ffopen(matfile, "w"); |
454 | write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2", |
455 | nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels); |
456 | gmx_ffclose(out); |
457 | } |
458 | |
459 | static void overlap(const char *outfile, int natoms, |
460 | rvec **eigvec1, |
461 | int nvec2, int *eignr2, rvec **eigvec2, |
462 | int noutvec, int *outvec, |
463 | const output_env_t oenv) |
464 | { |
465 | FILE *out; |
466 | int i, v, vec, x; |
467 | real overlap, inp; |
468 | |
469 | fprintf(stderrstderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n"); |
470 | for (i = 0; i < noutvec; i++) |
471 | { |
472 | fprintf(stderrstderr, "%d ", outvec[i]+1); |
473 | } |
474 | fprintf(stderrstderr, "\n"); |
475 | |
476 | out = xvgropen(outfile, "Subspace overlap", |
477 | "Eigenvectors of trajectory 2", "Overlap", oenv); |
478 | fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", |
479 | noutvec); |
480 | overlap = 0; |
481 | for (x = 0; x < nvec2; x++) |
482 | { |
483 | for (v = 0; v < noutvec; v++) |
484 | { |
485 | vec = outvec[v]; |
486 | inp = 0; |
487 | for (i = 0; i < natoms; i++) |
488 | { |
489 | inp += iprod(eigvec1[vec][i], eigvec2[x][i]); |
490 | } |
491 | overlap += sqr(inp); |
492 | } |
493 | fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec); |
494 | } |
495 | |
496 | gmx_ffclose(out); |
497 | } |
498 | |
499 | static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox, |
500 | const char *projfile, const char *twodplotfile, |
501 | const char *threedplotfile, const char *filterfile, int skip, |
502 | const char *extremefile, gmx_bool bExtrAll, real extreme, |
503 | int nextr, t_atoms *atoms, int natoms, atom_id *index, |
504 | gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls, |
505 | real *sqrtm, rvec *xav, |
506 | int *eignr, rvec **eigvec, |
507 | int noutvec, int *outvec, gmx_bool bSplit, |
508 | const output_env_t oenv) |
509 | { |
510 | FILE *xvgrout = NULL((void*)0); |
511 | int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame; |
512 | t_trxstatus *out = NULL((void*)0); |
513 | t_trxstatus *status; |
514 | int noutvec_extr, imin, imax; |
515 | real *pmin, *pmax; |
516 | atom_id *all_at; |
517 | matrix box; |
518 | rvec *xread, *x; |
519 | real t, inp, **inprod = NULL((void*)0), min = 0, max = 0; |
520 | char str[STRLEN4096], str2[STRLEN4096], **ylabel, *c; |
521 | real fact; |
522 | gmx_rmpbc_t gpbc = NULL((void*)0); |
523 | |
524 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 524, (natoms), sizeof(*(x))); |
525 | |
526 | if (bExtrAll) |
527 | { |
528 | noutvec_extr = noutvec; |
529 | } |
530 | else |
531 | { |
532 | noutvec_extr = 1; |
533 | } |
534 | |
535 | |
536 | if (trajfile) |
537 | { |
538 | snew(inprod, noutvec+1)(inprod) = save_calloc("inprod", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 538, (noutvec+1), sizeof(*(inprod))); |
539 | |
540 | if (filterfile) |
541 | { |
542 | fprintf(stderrstderr, "Writing a filtered trajectory to %s using eigenvectors\n", |
543 | filterfile); |
544 | for (i = 0; i < noutvec; i++) |
545 | { |
546 | fprintf(stderrstderr, "%d ", outvec[i]+1); |
547 | } |
548 | fprintf(stderrstderr, "\n"); |
549 | out = open_trx(filterfile, "w"); |
550 | } |
551 | snew_size = 0; |
552 | nfr = 0; |
553 | nframes = 0; |
554 | nat = read_first_x(oenv, &status, trajfile, &t, &xread, box); |
555 | if (nat > atoms->nr) |
556 | { |
557 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 557, "the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)", nat, atoms->nr); |
558 | } |
559 | snew(all_at, nat)(all_at) = save_calloc("all_at", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 559, (nat), sizeof(*(all_at))); |
560 | |
561 | if (top) |
562 | { |
563 | gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat); |
564 | } |
565 | |
566 | for (i = 0; i < nat; i++) |
567 | { |
568 | all_at[i] = i; |
569 | } |
570 | do |
571 | { |
572 | if (nfr % skip == 0) |
573 | { |
574 | if (top) |
575 | { |
576 | gmx_rmpbc(gpbc, nat, box, xread); |
577 | } |
578 | if (nframes >= snew_size) |
579 | { |
580 | snew_size += 100; |
581 | for (i = 0; i < noutvec+1; i++) |
582 | { |
583 | srenew(inprod[i], snew_size)(inprod[i]) = save_realloc("inprod[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 583, (inprod[i]), (snew_size), sizeof(*(inprod[i]))); |
584 | } |
585 | } |
586 | inprod[noutvec][nframes] = t; |
587 | /* calculate x: a fitted struture of the selected atoms */ |
588 | if (bFit) |
589 | { |
590 | reset_x(nfit, ifit, nat, NULL((void*)0), xread, w_rls); |
591 | do_fit(nat, w_rls, xref, xread); |
592 | } |
593 | for (i = 0; i < natoms; i++) |
594 | { |
595 | copy_rvec(xread[index[i]], x[i]); |
596 | } |
597 | |
598 | for (v = 0; v < noutvec; v++) |
599 | { |
600 | vec = outvec[v]; |
601 | /* calculate (mass-weighted) projection */ |
602 | inp = 0; |
603 | for (i = 0; i < natoms; i++) |
604 | { |
605 | inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+ |
606 | eigvec[vec][i][1]*(x[i][1]-xav[i][1])+ |
607 | eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i]; |
608 | } |
609 | inprod[v][nframes] = inp; |
610 | } |
611 | if (filterfile) |
612 | { |
613 | for (i = 0; i < natoms; i++) |
614 | { |
615 | for (d = 0; d < DIM3; d++) |
616 | { |
617 | /* misuse xread for output */ |
618 | xread[index[i]][d] = xav[i][d]; |
619 | for (v = 0; v < noutvec; v++) |
620 | { |
621 | xread[index[i]][d] += |
622 | inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i]; |
623 | } |
624 | } |
625 | } |
626 | write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL((void*)0), NULL((void*)0)); |
627 | } |
628 | nframes++; |
629 | } |
630 | nfr++; |
631 | } |
632 | while (read_next_x(oenv, status, &t, xread, box)); |
633 | close_trx(status); |
634 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 634, (x)); |
635 | if (filterfile) |
636 | { |
637 | close_trx(out); |
638 | } |
639 | } |
640 | else |
641 | { |
642 | snew(xread, atoms->nr)(xread) = save_calloc("xread", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 642, (atoms->nr), sizeof(*(xread))); |
643 | } |
644 | |
645 | if (top) |
646 | { |
647 | gmx_rmpbc_done(gpbc); |
648 | } |
649 | |
650 | |
651 | if (projfile) |
652 | { |
653 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 653, (noutvec), sizeof(*(ylabel))); |
654 | for (v = 0; v < noutvec; v++) |
655 | { |
656 | sprintf(str, "vec %d", eignr[outvec[v]]+1); |
657 | ylabel[v] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); |
658 | } |
659 | sprintf(str, "projection on eigenvectors (%s)", proj_unit); |
660 | write_xvgr_graphs(projfile, noutvec, 1, str, NULL((void*)0), output_env_get_xvgr_tlabel(oenv), |
661 | (const char **)ylabel, |
662 | nframes, inprod[noutvec], inprod, NULL((void*)0), |
663 | output_env_get_time_factor(oenv), FALSE0, bSplit, oenv); |
664 | } |
665 | |
666 | if (twodplotfile) |
667 | { |
668 | sprintf(str, "projection on eigenvector %d (%s)", |
669 | eignr[outvec[0]]+1, proj_unit); |
670 | sprintf(str2, "projection on eigenvector %d (%s)", |
671 | eignr[outvec[noutvec-1]]+1, proj_unit); |
672 | xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2, |
673 | oenv); |
674 | for (i = 0; i < nframes; i++) |
675 | { |
676 | if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5) |
677 | { |
678 | fprintf(xvgrout, "&\n"); |
679 | } |
680 | fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]); |
681 | } |
682 | gmx_ffclose(xvgrout); |
683 | } |
684 | |
685 | if (threedplotfile) |
686 | { |
687 | t_atoms atoms; |
688 | rvec *x; |
689 | real *b = NULL((void*)0); |
690 | matrix box; |
691 | char *resnm, *atnm, pdbform[STRLEN4096]; |
692 | gmx_bool bPDB, b4D; |
693 | FILE *out; |
694 | |
695 | if (noutvec < 3) |
696 | { |
697 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 697, "You have selected less than 3 eigenvectors"); |
698 | } |
699 | |
700 | /* initialize */ |
701 | bPDB = fn2ftp(threedplotfile) == efPDB; |
702 | clear_mat(box); |
703 | box[XX0][XX0] = box[YY1][YY1] = box[ZZ2][ZZ2] = 1; |
704 | |
705 | b4D = bPDB && (noutvec >= 4); |
706 | if (b4D) |
707 | { |
708 | fprintf(stderrstderr, "You have selected four or more eigenvectors:\n" |
709 | "fourth eigenvector will be plotted " |
710 | "in bfactor field of pdb file\n"); |
711 | sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d", |
712 | eignr[outvec[0]]+1, eignr[outvec[1]]+1, |
713 | eignr[outvec[2]]+1, eignr[outvec[3]]+1); |
714 | } |
715 | else |
716 | { |
717 | sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d", |
718 | eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1); |
719 | } |
720 | init_t_atoms(&atoms, nframes, FALSE0); |
721 | snew(x, nframes)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 721, (nframes), sizeof(*(x))); |
722 | snew(b, nframes)(b) = save_calloc("b", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 722, (nframes), sizeof(*(b))); |
723 | atnm = strdup("C")(__extension__ (__builtin_constant_p ("C") && ((size_t )(const void *)(("C") + 1) - (size_t)(const void *)("C") == 1 ) ? (((const char *) ("C"))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen ("C") + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, "C", __len); __retval ; })) : __strdup ("C"))); |
724 | resnm = strdup("PRJ")(__extension__ (__builtin_constant_p ("PRJ") && ((size_t )(const void *)(("PRJ") + 1) - (size_t)(const void *)("PRJ") == 1) ? (((const char *) ("PRJ"))[0] == '\0' ? (char *) calloc ( (size_t) 1, (size_t) 1) : ({ size_t __len = strlen ("PRJ") + 1 ; char *__retval = (char *) malloc (__len); if (__retval != ( (void*)0)) __retval = (char *) memcpy (__retval, "PRJ", __len ); __retval; })) : __strdup ("PRJ"))); |
725 | |
726 | if (nframes > 10000) |
727 | { |
728 | fact = 10000.0/nframes; |
729 | } |
730 | else |
731 | { |
732 | fact = 1.0; |
733 | } |
734 | |
735 | for (i = 0; i < nframes; i++) |
736 | { |
737 | atoms.atomname[i] = &atnm; |
738 | atoms.atom[i].resind = i; |
739 | atoms.resinfo[i].name = &resnm; |
740 | atoms.resinfo[i].nr = ceil(i*fact); |
741 | atoms.resinfo[i].ic = ' '; |
742 | x[i][XX0] = inprod[0][i]; |
743 | x[i][YY1] = inprod[1][i]; |
744 | x[i][ZZ2] = inprod[2][i]; |
745 | if (b4D) |
746 | { |
747 | b[i] = inprod[3][i]; |
748 | } |
749 | } |
750 | if ( ( b4D || bSplit ) && bPDB) |
751 | { |
752 | strcpy(pdbform, get_pdbformat()); |
753 | strcat(pdbform, "%8.4f%8.4f\n"); |
754 | |
755 | out = gmx_ffopen(threedplotfile, "w"); |
756 | fprintf(out, "HEADER %s\n", str); |
757 | if (b4D) |
758 | { |
759 | fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor"); |
760 | } |
761 | j = 0; |
762 | for (i = 0; i < atoms.nr; i++) |
763 | { |
764 | if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5) |
765 | { |
766 | fprintf(out, "TER\n"); |
767 | j = 0; |
768 | } |
769 | fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1, |
770 | 10*x[i][XX0], 10*x[i][YY1], 10*x[i][ZZ2], 1.0, 10*b[i]); |
771 | if (j > 0) |
772 | { |
773 | fprintf(out, "CONECT%5d%5d\n", i, i+1); |
774 | } |
775 | j++; |
776 | } |
777 | fprintf(out, "TER\n"); |
778 | gmx_ffclose(out); |
779 | } |
780 | else |
781 | { |
782 | write_sto_conf(threedplotfile, str, &atoms, x, NULL((void*)0), ePBC, box); |
783 | } |
784 | free_t_atoms(&atoms, FALSE0); |
785 | } |
786 | |
787 | if (extremefile) |
788 | { |
789 | snew(pmin, noutvec_extr)(pmin) = save_calloc("pmin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 789, (noutvec_extr), sizeof(*(pmin))); |
790 | snew(pmax, noutvec_extr)(pmax) = save_calloc("pmax", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 790, (noutvec_extr), sizeof(*(pmax))); |
791 | if (extreme == 0) |
792 | { |
793 | fprintf(stderrstderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum"); |
794 | fprintf(stderrstderr, |
795 | "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame"); |
796 | imin = 0; |
797 | imax = 0; |
798 | for (v = 0; v < noutvec_extr; v++) |
799 | { |
800 | for (i = 0; i < nframes; i++) |
801 | { |
802 | if (inprod[v][i] < inprod[v][imin]) |
803 | { |
804 | imin = i; |
805 | } |
806 | if (inprod[v][i] > inprod[v][imax]) |
807 | { |
808 | imax = i; |
809 | } |
810 | } |
811 | pmin[v] = inprod[v][imin]; |
812 | pmax[v] = inprod[v][imax]; |
813 | fprintf(stderrstderr, "%7d %10.6f %10d %10.6f %10d\n", |
814 | eignr[outvec[v]]+1, |
815 | pmin[v], imin, pmax[v], imax); |
816 | } |
817 | } |
818 | else |
819 | { |
820 | pmin[0] = -extreme; |
821 | pmax[0] = extreme; |
822 | } |
823 | /* build format string for filename: */ |
824 | strcpy(str, extremefile); /* copy filename */ |
825 | c = strrchr(str, '.'); /* find where extention begins */ |
826 | strcpy(str2, c); /* get extention */ |
827 | sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */ |
828 | for (v = 0; v < noutvec_extr; v++) |
829 | { |
830 | /* make filename using format string */ |
831 | if (noutvec_extr == 1) |
832 | { |
833 | strcpy(str2, extremefile); |
834 | } |
835 | else |
836 | { |
837 | sprintf(str2, str, eignr[outvec[v]]+1); |
838 | } |
839 | fprintf(stderrstderr, "Writing %d frames along eigenvector %d to %s\n", |
840 | nextr, outvec[v]+1, str2); |
841 | out = open_trx(str2, "w"); |
842 | for (frame = 0; frame < nextr; frame++) |
843 | { |
844 | if ((extreme == 0) && (nextr <= 3)) |
845 | { |
846 | for (i = 0; i < natoms; i++) |
847 | { |
848 | atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame; |
849 | } |
850 | } |
851 | for (i = 0; i < natoms; i++) |
852 | { |
853 | for (d = 0; d < DIM3; d++) |
854 | { |
855 | xread[index[i]][d] = |
856 | (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1) |
857 | *eigvec[outvec[v]][i][d]/sqrtm[i]); |
858 | } |
859 | } |
860 | write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL((void*)0), NULL((void*)0)); |
861 | } |
862 | close_trx(out); |
863 | } |
864 | sfree(pmin)save_free("pmin", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 864, (pmin)); |
865 | sfree(pmax)save_free("pmax", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 865, (pmax)); |
866 | } |
867 | fprintf(stderrstderr, "\n"); |
868 | } |
869 | |
870 | static void components(const char *outfile, int natoms, |
871 | int *eignr, rvec **eigvec, |
872 | int noutvec, int *outvec, |
873 | const output_env_t oenv) |
874 | { |
875 | int g, s, v, i; |
876 | real *x, ***y; |
877 | char str[STRLEN4096], **ylabel; |
878 | |
879 | fprintf(stderrstderr, "Writing eigenvector components to %s\n", outfile); |
880 | |
881 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 881, (noutvec), sizeof(*(ylabel))); |
882 | snew(y, noutvec)(y) = save_calloc("y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 882, (noutvec), sizeof(*(y))); |
883 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 883, (natoms), sizeof(*(x))); |
884 | for (i = 0; i < natoms; i++) |
885 | { |
886 | x[i] = i+1; |
887 | } |
888 | for (g = 0; g < noutvec; g++) |
889 | { |
890 | v = outvec[g]; |
891 | sprintf(str, "vec %d", eignr[v]+1); |
892 | ylabel[g] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); |
893 | snew(y[g], 4)(y[g]) = save_calloc("y[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 893, (4), sizeof(*(y[g]))); |
894 | for (s = 0; s < 4; s++) |
895 | { |
896 | snew(y[g][s], natoms)(y[g][s]) = save_calloc("y[g][s]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 896, (natoms), sizeof(*(y[g][s]))); |
897 | } |
898 | for (i = 0; i < natoms; i++) |
899 | { |
900 | y[g][0][i] = norm(eigvec[v][i]); |
901 | for (s = 0; s < 3; s++) |
902 | { |
903 | y[g][s+1][i] = eigvec[v][i][s]; |
904 | } |
905 | } |
906 | } |
907 | write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components", |
908 | "black: total, red: x, green: y, blue: z", |
909 | "Atom number", (const char **)ylabel, |
910 | natoms, x, NULL((void*)0), y, 1, FALSE0, FALSE0, oenv); |
911 | fprintf(stderrstderr, "\n"); |
912 | } |
913 | |
914 | static void rmsf(const char *outfile, int natoms, real *sqrtm, |
915 | int *eignr, rvec **eigvec, |
916 | int noutvec, int *outvec, |
917 | real *eigval, int neig, |
918 | const output_env_t oenv) |
919 | { |
920 | int nrow, g, v, i; |
921 | real *x, **y; |
922 | char str[STRLEN4096], **ylabel; |
923 | |
924 | for (i = 0; i < neig; i++) |
925 | { |
926 | if (eigval[i] < 0) |
927 | { |
928 | eigval[i] = 0; |
929 | } |
930 | } |
931 | |
932 | fprintf(stderrstderr, "Writing rmsf to %s\n", outfile); |
933 | |
934 | snew(ylabel, noutvec)(ylabel) = save_calloc("ylabel", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 934, (noutvec), sizeof(*(ylabel))); |
935 | snew(y, noutvec)(y) = save_calloc("y", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 935, (noutvec), sizeof(*(y))); |
936 | snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 936, (natoms), sizeof(*(x))); |
937 | for (i = 0; i < natoms; i++) |
938 | { |
939 | x[i] = i+1; |
940 | } |
941 | for (g = 0; g < noutvec; g++) |
942 | { |
943 | v = outvec[g]; |
944 | if (eignr[v] >= neig) |
945 | { |
946 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 946, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig); |
947 | } |
948 | sprintf(str, "vec %d", eignr[v]+1); |
949 | ylabel[g] = strdup(str)(__extension__ (__builtin_constant_p (str) && ((size_t )(const void *)((str) + 1) - (size_t)(const void *)(str) == 1 ) ? (((const char *) (str))[0] == '\0' ? (char *) calloc ((size_t ) 1, (size_t) 1) : ({ size_t __len = strlen (str) + 1; char * __retval = (char *) malloc (__len); if (__retval != ((void*)0 )) __retval = (char *) memcpy (__retval, str, __len); __retval ; })) : __strdup (str))); |
950 | snew(y[g], natoms)(y[g]) = save_calloc("y[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 950, (natoms), sizeof(*(y[g]))); |
951 | for (i = 0; i < natoms; i++) |
952 | { |
953 | y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i]; |
954 | } |
955 | } |
956 | write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL((void*)0), |
957 | "Atom number", (const char **)ylabel, |
958 | natoms, x, y, NULL((void*)0), 1, TRUE1, FALSE0, oenv); |
959 | fprintf(stderrstderr, "\n"); |
960 | } |
961 | |
962 | int gmx_anaeig(int argc, char *argv[]) |
963 | { |
964 | static const char *desc[] = { |
965 | "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a", |
966 | "covariance matrix ([gmx-covar]) or of a Normal Modes analysis", |
967 | "([gmx-nmeig]).[PAR]", |
968 | |
969 | "When a trajectory is projected on eigenvectors, all structures are", |
970 | "fitted to the structure in the eigenvector file, if present, otherwise", |
971 | "to the structure in the structure file. When no run input file is", |
972 | "supplied, periodicity will not be taken into account. Most analyses", |
973 | "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when", |
974 | "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]", |
975 | |
976 | "[TT]-comp[tt]: plot the vector components per atom of eigenvectors", |
977 | "[TT]-first[tt] to [TT]-last[tt].[PAR]", |
978 | |
979 | "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors", |
980 | "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]", |
981 | |
982 | "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors", |
983 | "[TT]-first[tt] to [TT]-last[tt].", |
984 | "The projections of a trajectory on the eigenvectors of its", |
985 | "covariance matrix are called principal components (pc's).", |
986 | "It is often useful to check the cosine content of the pc's,", |
987 | "since the pc's of random diffusion are cosines with the number", |
988 | "of periods equal to half the pc index.", |
989 | "The cosine content of the pc's can be calculated with the program", |
990 | "[gmx-analyze].[PAR]", |
991 | |
992 | "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors", |
993 | "[TT]-first[tt] and [TT]-last[tt].[PAR]", |
994 | |
995 | "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first", |
996 | "three selected eigenvectors.[PAR]", |
997 | |
998 | "[TT]-filt[tt]: filter the trajectory to show only the motion along", |
999 | "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]", |
1000 | |
1001 | "[TT]-extr[tt]: calculate the two extreme projections along a trajectory", |
1002 | "on the average structure and interpolate [TT]-nframes[tt] frames", |
1003 | "between them, or set your own extremes with [TT]-max[tt]. The", |
1004 | "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and", |
1005 | "[TT]-last[tt] have been set explicitly, in which case all eigenvectors", |
1006 | "will be written to separate files. Chain identifiers will be added", |
1007 | "when writing a [TT].pdb[tt] file with two or three structures (you", |
1008 | "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]", |
1009 | |
1010 | " Overlap calculations between covariance analysis:[BR]", |
1011 | " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]", |
1012 | |
1013 | "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in", |
1014 | "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]", |
1015 | "in file [TT]-v[tt].[PAR]", |
1016 | |
1017 | "[TT]-inpr[tt]: calculate a matrix of inner-products between", |
1018 | "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors", |
1019 | "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]", |
1020 | "have been set explicitly.[PAR]", |
1021 | |
1022 | "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,", |
1023 | "a single number for the overlap between the covariance matrices is", |
1024 | "generated. The formulas are:[BR]", |
1025 | " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]", |
1026 | "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]", |
1027 | " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]", |
1028 | "where M1 and M2 are the two covariance matrices and tr is the trace", |
1029 | "of a matrix. The numbers are proportional to the overlap of the square", |
1030 | "root of the fluctuations. The normalized overlap is the most useful", |
1031 | "number, it is 1 for identical matrices and 0 when the sampled", |
1032 | "subspaces are orthogonal.[PAR]", |
1033 | "When the [TT]-entropy[tt] flag is given an entropy estimate will be", |
1034 | "computed based on the Quasiharmonic approach and based on", |
1035 | "Schlitter's formula." |
1036 | }; |
1037 | static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6; |
1038 | static real max = 0.0, temp = 298.15; |
1039 | static gmx_bool bSplit = FALSE0, bEntropy = FALSE0; |
1040 | t_pargs pa[] = { |
1041 | { "-first", FALSE0, etINT, {&first}, |
1042 | "First eigenvector for analysis (-1 is select)" }, |
1043 | { "-last", FALSE0, etINT, {&last}, |
1044 | "Last eigenvector for analysis (-1 is till the last)" }, |
1045 | { "-skip", FALSE0, etINT, {&skip}, |
1046 | "Only analyse every nr-th frame" }, |
1047 | { "-max", FALSE0, etREAL, {&max}, |
1048 | "Maximum for projection of the eigenvector on the average structure, " |
1049 | "max=0 gives the extremes" }, |
1050 | { "-nframes", FALSE0, etINT, {&nextr}, |
1051 | "Number of frames for the extremes output" }, |
1052 | { "-split", FALSE0, etBOOL, {&bSplit}, |
1053 | "Split eigenvector projections where time is zero" }, |
1054 | { "-entropy", FALSE0, etBOOL, {&bEntropy}, |
1055 | "Compute entropy according to the Quasiharmonic formula or Schlitter's method." }, |
1056 | { "-temp", FALSE0, etREAL, {&temp}, |
1057 | "Temperature for entropy calculations" }, |
1058 | { "-nevskip", FALSE0, etINT, {&nskip}, |
1059 | "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." } |
1060 | }; |
1061 | #define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))) |
1062 | |
1063 | FILE *out; |
1064 | int status, trjout; |
1065 | t_topology top; |
1066 | int ePBC = -1; |
1067 | t_atoms *atoms = NULL((void*)0); |
1068 | rvec *xtop, *xref1, *xref2, *xrefp = NULL((void*)0); |
1069 | gmx_bool bDMR1, bDMA1, bDMR2, bDMA2; |
1070 | int nvec1, nvec2, *eignr1 = NULL((void*)0), *eignr2 = NULL((void*)0); |
1071 | rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL((void*)0), **eigvec2 = NULL((void*)0); |
1072 | matrix topbox; |
1073 | real xid, totmass, *sqrtm, *w_rls, t, lambda; |
1074 | int natoms, step; |
1075 | char *grpname; |
1076 | const char *indexfile; |
1077 | char title[STRLEN4096]; |
1078 | int i, j, d; |
1079 | int nout, *iout, noutvec, *outvec, nfit; |
1080 | atom_id *index, *ifit; |
1081 | const char *VecFile, *Vec2File, *topfile; |
1082 | const char *EigFile, *Eig2File; |
1083 | const char *CompFile, *RmsfFile, *ProjOnVecFile; |
1084 | const char *TwoDPlotFile, *ThreeDPlotFile; |
1085 | const char *FilterFile, *ExtremeFile; |
1086 | const char *OverlapFile, *InpMatFile; |
1087 | gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj; |
1088 | gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D; |
1089 | real *eigval1 = NULL((void*)0), *eigval2 = NULL((void*)0); |
1090 | int neig1, neig2; |
1091 | double **xvgdata; |
1092 | output_env_t oenv; |
1093 | gmx_rmpbc_t gpbc; |
1094 | |
1095 | t_filenm fnm[] = { |
1096 | { efTRN, "-v", "eigenvec", ffREAD1<<1 }, |
1097 | { efTRN, "-v2", "eigenvec2", ffOPTRD(1<<1 | 1<<3) }, |
1098 | { efTRX, "-f", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1099 | { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1100 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
1101 | { efXVG, "-eig", "eigenval", ffOPTRD(1<<1 | 1<<3) }, |
1102 | { efXVG, "-eig2", "eigenval2", ffOPTRD(1<<1 | 1<<3) }, |
1103 | { efXVG, "-comp", "eigcomp", ffOPTWR(1<<2| 1<<3) }, |
1104 | { efXVG, "-rmsf", "eigrmsf", ffOPTWR(1<<2| 1<<3) }, |
1105 | { efXVG, "-proj", "proj", ffOPTWR(1<<2| 1<<3) }, |
1106 | { efXVG, "-2d", "2dproj", ffOPTWR(1<<2| 1<<3) }, |
1107 | { efSTO, "-3d", "3dproj.pdb", ffOPTWR(1<<2| 1<<3) }, |
1108 | { efTRX, "-filt", "filtered", ffOPTWR(1<<2| 1<<3) }, |
1109 | { efTRX, "-extr", "extreme.pdb", ffOPTWR(1<<2| 1<<3) }, |
1110 | { efXVG, "-over", "overlap", ffOPTWR(1<<2| 1<<3) }, |
1111 | { efXPM, "-inpr", "inprod", ffOPTWR(1<<2| 1<<3) } |
1112 | }; |
1113 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
1114 | |
1115 | if (!parse_common_args(&argc, argv, |
1116 | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13), |
1117 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
1118 | { |
1119 | return 0; |
1120 | } |
1121 | |
1122 | indexfile = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1123 | |
1124 | VecFile = opt2fn("-v", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1125 | Vec2File = opt2fn_null("-v2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1126 | topfile = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1127 | EigFile = opt2fn_null("-eig", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1128 | Eig2File = opt2fn_null("-eig2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1129 | CompFile = opt2fn_null("-comp", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1130 | RmsfFile = opt2fn_null("-rmsf", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1131 | ProjOnVecFile = opt2fn_null("-proj", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1132 | TwoDPlotFile = opt2fn_null("-2d", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1133 | ThreeDPlotFile = opt2fn_null("-3d", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1134 | FilterFile = opt2fn_null("-filt", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1135 | ExtremeFile = opt2fn_null("-extr", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1136 | OverlapFile = opt2fn_null("-over", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1137 | InpMatFile = ftp2fn_null(efXPM, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1138 | |
1139 | bTop = fn2bTPX(topfile); |
1140 | bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile |
1141 | || FilterFile || ExtremeFile; |
1142 | bFirstLastSet = |
1143 | opt2parg_bSet("-first", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa) && opt2parg_bSet("-last", NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa); |
1144 | bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile || |
1145 | OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet); |
1146 | bVec2 = Vec2File || OverlapFile || InpMatFile; |
1147 | bM = RmsfFile || bProj; |
1148 | bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0)) |
1149 | || TwoDPlotFile || ThreeDPlotFile; |
1150 | bIndex = bM || bProj; |
1151 | bTPS = ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bM || bTraj || |
1152 | FilterFile || (bIndex && indexfile); |
1153 | bCompare = Vec2File || Eig2File; |
1154 | bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB; |
1155 | |
1156 | read_eigenvectors(VecFile, &natoms, &bFit1, |
1157 | &xref1, &bDMR1, &xav1, &bDMA1, |
1158 | &nvec1, &eignr1, &eigvec1, &eigval1); |
1159 | neig1 = DIM3*natoms; |
1160 | |
1161 | /* Overwrite eigenvalues from separate files if the user provides them */ |
1162 | if (EigFile != NULL((void*)0)) |
1163 | { |
1164 | int neig_tmp = read_xvg(EigFile, &xvgdata, &i); |
1165 | if (neig_tmp != neig1) |
1166 | { |
1167 | fprintf(stderrstderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms); |
1168 | } |
1169 | neig1 = neig_tmp; |
1170 | srenew(eigval1, neig1)(eigval1) = save_realloc("eigval1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1170, (eigval1), (neig1), sizeof(*(eigval1))); |
1171 | for (j = 0; j < neig1; j++) |
1172 | { |
1173 | real tmp = eigval1[j]; |
1174 | eigval1[j] = xvgdata[1][j]; |
1175 | if (debug && (eigval1[j] != tmp)) |
1176 | { |
1177 | fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n", |
1178 | j, tmp, eigval1[j]); |
1179 | } |
1180 | } |
1181 | for (j = 0; j < i; j++) |
1182 | { |
1183 | sfree(xvgdata[j])save_free("xvgdata[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1183, (xvgdata[j])); |
1184 | } |
1185 | sfree(xvgdata)save_free("xvgdata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1185, (xvgdata)); |
1186 | fprintf(stderrstderr, "Read %d eigenvalues from %s\n", neig1, EigFile); |
1187 | } |
1188 | |
1189 | if (bEntropy) |
1190 | { |
1191 | if (bDMA1) |
1192 | { |
1193 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1193, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting"); |
1194 | } |
1195 | calc_entropy_qh(stdoutstdout, neig1, eigval1, temp, nskip); |
1196 | calc_entropy_schlitter(stdoutstdout, neig1, nskip, eigval1, temp); |
1197 | } |
1198 | |
1199 | if (bVec2) |
1200 | { |
1201 | if (!Vec2File) |
1202 | { |
1203 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1203, "Need a second eigenvector file to do this analysis."); |
1204 | } |
1205 | read_eigenvectors(Vec2File, &neig2, &bFit2, |
1206 | &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2); |
1207 | |
1208 | neig2 = DIM3*neig2; |
1209 | if (neig2 != neig1) |
1210 | { |
1211 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1211, "Dimensions in the eigenvector files don't match"); |
1212 | } |
1213 | } |
1214 | |
1215 | if (Eig2File != NULL((void*)0)) |
1216 | { |
1217 | neig2 = read_xvg(Eig2File, &xvgdata, &i); |
1218 | srenew(eigval2, neig2)(eigval2) = save_realloc("eigval2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1218, (eigval2), (neig2), sizeof(*(eigval2))); |
1219 | for (j = 0; j < neig2; j++) |
1220 | { |
1221 | eigval2[j] = xvgdata[1][j]; |
1222 | } |
1223 | for (j = 0; j < i; j++) |
1224 | { |
1225 | sfree(xvgdata[j])save_free("xvgdata[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1225, (xvgdata[j])); |
1226 | } |
1227 | sfree(xvgdata)save_free("xvgdata", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1227, (xvgdata)); |
1228 | fprintf(stderrstderr, "Read %d eigenvalues from %s\n", neig2, Eig2File); |
1229 | } |
1230 | |
1231 | |
1232 | if ((!bFit1 || xref1) && !bDMR1 && !bDMA1) |
1233 | { |
1234 | bM = FALSE0; |
1235 | } |
1236 | if ((xref1 == NULL((void*)0)) && (bM || bTraj)) |
1237 | { |
1238 | bTPS = TRUE1; |
1239 | } |
1240 | |
1241 | xtop = NULL((void*)0); |
1242 | nfit = 0; |
1243 | ifit = NULL((void*)0); |
1244 | w_rls = NULL((void*)0); |
1245 | |
1246 | if (!bTPS) |
1247 | { |
1248 | bTop = FALSE0; |
1249 | } |
1250 | else |
1251 | { |
1252 | bTop = read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
1253 | title, &top, &ePBC, &xtop, NULL((void*)0), topbox, bM); |
1254 | atoms = &top.atoms; |
1255 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr); |
1256 | gmx_rmpbc(gpbc, atoms->nr, topbox, xtop); |
1257 | /* Fitting is only required for the projection */ |
1258 | if (bProj && bFit1) |
1259 | { |
1260 | if (xref1 == NULL((void*)0)) |
1261 | { |
1262 | printf("\nNote: the structure in %s should be the same\n" |
1263 | " as the one used for the fit in g_covar\n", topfile); |
1264 | } |
1265 | printf("\nSelect the index group that was used for the least squares fit in g_covar\n"); |
1266 | get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname); |
1267 | |
1268 | snew(w_rls, atoms->nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1268, (atoms->nr), sizeof(*(w_rls))); |
1269 | for (i = 0; (i < nfit); i++) |
1270 | { |
1271 | if (bDMR1) |
1272 | { |
1273 | w_rls[ifit[i]] = atoms->atom[ifit[i]].m; |
1274 | } |
1275 | else |
1276 | { |
1277 | w_rls[ifit[i]] = 1.0; |
1278 | } |
1279 | } |
1280 | |
1281 | snew(xrefp, atoms->nr)(xrefp) = save_calloc("xrefp", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1281, (atoms->nr), sizeof(*(xrefp))); |
1282 | if (xref1 != NULL((void*)0)) |
1283 | { |
1284 | /* Safety check between selected fit-group and reference structure read from the eigenvector file */ |
1285 | if (natoms != nfit) |
1286 | { |
1287 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1287, "you selected a group with %d elements instead of %d, your selection does not fit the reference structure in the eigenvector file.", nfit, natoms); |
1288 | } |
1289 | for (i = 0; (i < nfit); i++) |
1290 | { |
1291 | copy_rvec(xref1[i], xrefp[ifit[i]]); |
1292 | } |
1293 | } |
1294 | else |
1295 | { |
1296 | /* The top coordinates are the fitting reference */ |
1297 | for (i = 0; (i < nfit); i++) |
1298 | { |
1299 | copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]); |
1300 | } |
1301 | reset_x(nfit, ifit, atoms->nr, NULL((void*)0), xrefp, w_rls); |
1302 | } |
1303 | } |
1304 | gmx_rmpbc_done(gpbc); |
1305 | } |
1306 | |
1307 | if (bIndex) |
1308 | { |
1309 | printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms); |
1310 | get_index(atoms, indexfile, 1, &i, &index, &grpname); |
1311 | if (i != natoms) |
1312 | { |
1313 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1313, "you selected a group with %d elements instead of %d", i, natoms); |
1314 | } |
1315 | printf("\n"); |
1316 | } |
1317 | |
1318 | snew(sqrtm, natoms)(sqrtm) = save_calloc("sqrtm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1318, (natoms), sizeof(*(sqrtm))); |
1319 | if (bM && bDMA1) |
1320 | { |
1321 | proj_unit = "u\\S1/2\\Nnm"; |
1322 | for (i = 0; (i < natoms); i++) |
1323 | { |
1324 | sqrtm[i] = sqrt(atoms->atom[index[i]].m); |
1325 | } |
1326 | } |
1327 | else |
1328 | { |
1329 | proj_unit = "nm"; |
1330 | for (i = 0; (i < natoms); i++) |
1331 | { |
1332 | sqrtm[i] = 1.0; |
1333 | } |
1334 | } |
1335 | |
1336 | if (bVec2) |
1337 | { |
1338 | t = 0; |
1339 | totmass = 0; |
1340 | for (i = 0; (i < natoms); i++) |
1341 | { |
1342 | for (d = 0; (d < DIM3); d++) |
1343 | { |
1344 | t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]); |
1345 | totmass += sqr(sqrtm[i]); |
1346 | } |
1347 | } |
1348 | fprintf(stdoutstdout, "RMSD (without fit) between the two average structures:" |
1349 | " %.3f (nm)\n\n", sqrt(t/totmass)); |
1350 | } |
1351 | |
1352 | if (last == -1) |
1353 | { |
1354 | last = natoms*DIM3; |
1355 | } |
1356 | if (first > -1) |
1357 | { |
1358 | if (bFirstToLast) |
1359 | { |
1360 | /* make an index from first to last */ |
1361 | nout = last-first+1; |
1362 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1362, (nout), sizeof(*(iout))); |
1363 | for (i = 0; i < nout; i++) |
1364 | { |
1365 | iout[i] = first-1+i; |
1366 | } |
1367 | } |
1368 | else if (ThreeDPlotFile) |
1369 | { |
1370 | /* make an index of first+(0,1,2) and last */ |
1371 | nout = bPDB3D ? 4 : 3; |
1372 | nout = min(last-first+1, nout)(((last-first+1) < (nout)) ? (last-first+1) : (nout) ); |
1373 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1373, (nout), sizeof(*(iout))); |
1374 | iout[0] = first-1; |
1375 | iout[1] = first; |
1376 | if (nout > 3) |
1377 | { |
1378 | iout[2] = first+1; |
1379 | } |
1380 | iout[nout-1] = last-1; |
1381 | } |
1382 | else |
1383 | { |
1384 | /* make an index of first and last */ |
1385 | nout = 2; |
1386 | snew(iout, nout)(iout) = save_calloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1386, (nout), sizeof(*(iout))); |
1387 | iout[0] = first-1; |
1388 | iout[1] = last-1; |
1389 | } |
1390 | } |
1391 | else |
1392 | { |
1393 | printf("Select eigenvectors for output, end your selection with 0\n"); |
1394 | nout = -1; |
1395 | iout = NULL((void*)0); |
1396 | |
1397 | do |
1398 | { |
1399 | nout++; |
1400 | srenew(iout, nout+1)(iout) = save_realloc("iout", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1400, (iout), (nout+1), sizeof(*(iout))); |
1401 | if (1 != scanf("%d", &iout[nout])) |
1402 | { |
1403 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1403, "Error reading user input"); |
1404 | } |
1405 | iout[nout]--; |
1406 | } |
1407 | while (iout[nout] >= 0); |
1408 | |
1409 | printf("\n"); |
1410 | } |
1411 | /* make an index of the eigenvectors which are present */ |
1412 | snew(outvec, nout)(outvec) = save_calloc("outvec", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anaeig.c" , 1412, (nout), sizeof(*(outvec))); |
1413 | noutvec = 0; |
1414 | for (i = 0; i < nout; i++) |
1415 | { |
1416 | j = 0; |
1417 | while ((j < nvec1) && (eignr1[j] != iout[i])) |
1418 | { |
1419 | j++; |
1420 | } |
1421 | if ((j < nvec1) && (eignr1[j] == iout[i])) |
1422 | { |
1423 | outvec[noutvec] = j; |
1424 | noutvec++; |
1425 | } |
1426 | } |
1427 | fprintf(stderrstderr, "%d eigenvectors selected for output", noutvec); |
1428 | if (noutvec <= 100) |
1429 | { |
1430 | fprintf(stderrstderr, ":"); |
1431 | for (j = 0; j < noutvec; j++) |
1432 | { |
1433 | fprintf(stderrstderr, " %d", eignr1[outvec[j]]+1); |
1434 | } |
1435 | } |
1436 | fprintf(stderrstderr, "\n"); |
1437 | |
1438 | if (CompFile) |
1439 | { |
1440 | components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv); |
1441 | } |
1442 | |
1443 | if (RmsfFile) |
1444 | { |
1445 | rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1, |
1446 | neig1, oenv); |
1447 | } |
1448 | |
1449 | if (bProj) |
1450 | { |
1451 | project(bTraj ? opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) : NULL((void*)0), |
1452 | bTop ? &top : NULL((void*)0), ePBC, topbox, |
1453 | ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip, |
1454 | ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index, |
1455 | bFit1, xrefp, nfit, ifit, w_rls, |
1456 | sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit, |
1457 | oenv); |
1458 | } |
1459 | |
1460 | if (OverlapFile) |
1461 | { |
1462 | overlap(OverlapFile, natoms, |
1463 | eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv); |
1464 | } |
1465 | |
1466 | if (InpMatFile) |
1467 | { |
1468 | inprod_matrix(InpMatFile, natoms, |
1469 | nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2, |
1470 | bFirstLastSet, noutvec, outvec); |
1471 | } |
1472 | |
1473 | if (bCompare) |
1474 | { |
1475 | compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2); |
1476 | } |
1477 | |
1478 | |
1479 | if (!CompFile && !bProj && !OverlapFile && !InpMatFile && |
1480 | !bCompare && !bEntropy) |
1481 | { |
1482 | fprintf(stderrstderr, "\nIf you want some output," |
1483 | " set one (or two or ...) of the output file options\n"); |
1484 | } |
1485 | |
1486 | |
1487 | view_all(oenv, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1488 | |
1489 | return 0; |
1490 | } |