File: | gromacs/gmxana/gmx_anaeig.c |
Location: | line 1475, column 9 |
Description: | Function call argument is an uninitialized value |
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; | |||
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 | } |