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