2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
45 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/futil.h"
53 #include "gromacs/fileio/pdbio.h"
54 #include "gromacs/fileio/confio.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/matio.h"
58 #include "gromacs/fileio/xvgr.h"
60 #include "gromacs/pbcutil/rmpbc.h"
63 #include "gromacs/math/units.h"
66 #include "gromacs/math/do_fit.h"
68 static void calc_entropy_qh(FILE *fp
, int n
, real eigval
[], real temp
, int nskip
)
71 double hwkT
, w
, dS
, S
= 0;
74 hbar
= PLANCK1
/(2*M_PI
);
75 for (i
= 0; (i
< n
-nskip
); i
++)
79 lambda
= eigval
[i
]*AMU
;
80 w
= sqrt(BOLTZMANN
*temp
/lambda
)/NANO
;
81 hwkT
= (hbar
*w
)/(BOLTZMANN
*temp
);
82 dS
= (hwkT
/(exp(hwkT
) - 1) - log(1-exp(-hwkT
)));
86 fprintf(debug
, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
87 i
, w
, lambda
, hwkT
, dS
);
92 fprintf(stderr
, "eigval[%d] = %g\n", i
, eigval
[i
]);
96 fprintf(fp
, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
100 static void calc_entropy_schlitter(FILE *fp
, int n
, int nskip
,
101 real
*eigval
, real temp
)
107 double hbar
, kt
, kteh
, S
;
109 hbar
= PLANCK1
/(2*M_PI
);
111 kteh
= kt
*exp(2.0)/(hbar
*hbar
)*AMU
*(NANO
*NANO
);
114 fprintf(debug
, "n = %d, nskip = %d kteh = %g\n", n
, nskip
, kteh
);
118 for (i
= 0; (i
< n
-nskip
); i
++)
120 dd
= 1+kteh
*eigval
[i
];
125 fprintf(fp
, "The Entropy due to the Schlitter formula is %g J/mol K\n", S
);
128 const char *proj_unit
;
130 static real
tick_spacing(real range
, int minticks
)
139 sp
= 0.2*exp(log(10)*ceil(log(range
)/log(10)));
140 while (range
/sp
< minticks
-1)
148 static void write_xvgr_graphs(const char *file
, int ngraphs
, int nsetspergraph
,
149 const char *title
, const char *subtitle
,
150 const char *xlabel
, const char **ylabel
,
151 int n
, real
*x
, real
**y
, real
***sy
,
152 real scale_x
, gmx_bool bZero
, gmx_bool bSplit
,
153 const output_env_t oenv
)
157 real min
, max
, xsp
, ysp
;
159 out
= gmx_ffopen(file
, "w");
160 if (output_env_get_xvg_format(oenv
) == exvgXMGRACE
)
162 fprintf(out
, "@ autoscale onread none\n");
164 for (g
= 0; g
< ngraphs
; g
++)
170 for (i
= 0; i
< n
; i
++)
186 for (s
= 0; s
< nsetspergraph
; s
++)
188 for (i
= 0; i
< n
; i
++)
190 if (sy
[g
][s
][i
] < min
)
194 if (sy
[g
][s
][i
] > max
)
207 min
= min
-0.1*(max
-min
);
209 max
= max
+0.1*(max
-min
);
210 xsp
= tick_spacing((x
[n
-1]-x
[0])*scale_x
, 4);
211 ysp
= tick_spacing(max
-min
, 3);
212 if (output_env_get_print_xvgr_codes(oenv
))
214 fprintf(out
, "@ with g%d\n@ g%d on\n", g
, g
);
217 fprintf(out
, "@ title \"%s\"\n", title
);
220 fprintf(out
, "@ subtitle \"%s\"\n", subtitle
);
225 fprintf(out
, "@ xaxis label \"%s\"\n", xlabel
);
229 fprintf(out
, "@ xaxis ticklabel off\n");
233 fprintf(out
, "@ world xmin %g\n", x
[0]*scale_x
);
234 fprintf(out
, "@ world xmax %g\n", x
[n
-1]*scale_x
);
235 fprintf(out
, "@ world ymin %g\n", min
);
236 fprintf(out
, "@ world ymax %g\n", max
);
238 fprintf(out
, "@ view xmin 0.15\n");
239 fprintf(out
, "@ view xmax 0.85\n");
240 fprintf(out
, "@ view ymin %g\n", 0.15+(ngraphs
-1-g
)*0.7/ngraphs
);
241 fprintf(out
, "@ view ymax %g\n", 0.15+(ngraphs
-g
)*0.7/ngraphs
);
242 fprintf(out
, "@ yaxis label \"%s\"\n", ylabel
[g
]);
243 fprintf(out
, "@ xaxis tick major %g\n", xsp
);
244 fprintf(out
, "@ xaxis tick minor %g\n", xsp
/2);
245 fprintf(out
, "@ xaxis ticklabel start type spec\n");
246 fprintf(out
, "@ xaxis ticklabel start %g\n", ceil(min
/xsp
)*xsp
);
247 fprintf(out
, "@ yaxis tick major %g\n", ysp
);
248 fprintf(out
, "@ yaxis tick minor %g\n", ysp
/2);
249 fprintf(out
, "@ yaxis ticklabel start type spec\n");
250 fprintf(out
, "@ yaxis ticklabel start %g\n", ceil(min
/ysp
)*ysp
);
251 if ((min
< 0) && (max
> 0))
253 fprintf(out
, "@ zeroxaxis bar on\n");
254 fprintf(out
, "@ zeroxaxis bar linestyle 3\n");
257 for (s
= 0; s
< nsetspergraph
; s
++)
259 for (i
= 0; i
< n
; i
++)
261 if (bSplit
&& i
> 0 && fabs(x
[i
]) < 1e-5)
263 if (output_env_get_print_xvgr_codes(oenv
))
268 fprintf(out
, "%10.4f %10.5f\n",
269 x
[i
]*scale_x
, y
? y
[g
][i
] : sy
[g
][s
][i
]);
271 if (output_env_get_print_xvgr_codes(oenv
))
281 compare(int natoms
, int n1
, rvec
**eigvec1
, int n2
, rvec
**eigvec2
,
282 real
*eigval1
, int neig1
, real
*eigval2
, int neig2
)
286 double sum1
, sum2
, trace1
, trace2
, sab
, samsb2
, tmp
, ip
;
290 n
= min(n
, min(neig1
, neig2
));
291 fprintf(stdout
, "Will compare the covariance matrices using %d dimensions\n", n
);
294 for (i
= 0; i
< n
; i
++)
301 eigval1
[i
] = sqrt(eigval1
[i
]);
304 for (i
= n
; i
< neig1
; i
++)
306 trace1
+= eigval1
[i
];
309 for (i
= 0; i
< n
; i
++)
316 eigval2
[i
] = sqrt(eigval2
[i
]);
319 for (i
= n
; i
< neig2
; i
++)
321 trace2
+= eigval2
[i
];
324 fprintf(stdout
, "Trace of the two matrices: %g and %g\n", sum1
, sum2
);
325 if (neig1
!= n
|| neig2
!= n
)
327 fprintf(stdout
, "this is %d%% and %d%% of the total trace\n",
328 (int)(100*sum1
/trace1
+0.5), (int)(100*sum2
/trace2
+0.5));
330 fprintf(stdout
, "Square root of the traces: %g and %g\n",
331 sqrt(sum1
), sqrt(sum2
));
334 for (i
= 0; i
< n
; i
++)
337 for (j
= 0; j
< n
; j
++)
340 for (k
= 0; k
< natoms
; k
++)
342 ip
+= iprod(eigvec1
[i
][k
], eigvec2
[j
][k
]);
344 tmp
+= eigval2
[j
]*ip
*ip
;
346 sab
+= eigval1
[i
]*tmp
;
349 samsb2
= sum1
+sum2
-2*sab
;
355 fprintf(stdout
, "The overlap of the covariance matrices:\n");
356 fprintf(stdout
, " normalized: %.3f\n", 1-sqrt(samsb2
/(sum1
+sum2
)));
357 tmp
= 1-sab
/sqrt(sum1
*sum2
);
362 fprintf(stdout
, " shape: %.3f\n", 1-sqrt(tmp
));
366 static void inprod_matrix(const char *matfile
, int natoms
,
367 int nvec1
, int *eignr1
, rvec
**eigvec1
,
368 int nvec2
, int *eignr2
, rvec
**eigvec2
,
369 gmx_bool bSelect
, int noutvec
, int *outvec
)
373 int i
, x1
, y1
, x
, y
, nlevels
;
375 real inp
, *t_x
, *t_y
, max
;
383 for (y1
= 0; y1
< nx
; y1
++)
385 if (outvec
[y1
] < nvec2
)
387 t_y
[ny
] = eignr2
[outvec
[y1
]]+1;
396 for (y
= 0; y
< ny
; y
++)
398 t_y
[y
] = eignr2
[y
]+1;
402 fprintf(stderr
, "Calculating inner-product matrix of %dx%d eigenvectors\n",
408 for (x1
= 0; x1
< nx
; x1
++)
419 t_x
[x1
] = eignr1
[x
]+1;
420 fprintf(stderr
, " %d", eignr1
[x
]+1);
421 for (y1
= 0; y1
< ny
; y1
++)
426 while (outvec
[y1
] >= nvec2
)
436 for (i
= 0; i
< natoms
; i
++)
438 inp
+= iprod(eigvec1
[x
][i
], eigvec2
[y
][i
]);
440 mat
[x1
][y1
] = fabs(inp
);
441 if (mat
[x1
][y1
] > max
)
447 fprintf(stderr
, "\n");
448 rlo
.r
= 1; rlo
.g
= 1; rlo
.b
= 1;
449 rhi
.r
= 0; rhi
.g
= 0; rhi
.b
= 0;
451 out
= gmx_ffopen(matfile
, "w");
452 write_xpm(out
, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
453 nx
, ny
, t_x
, t_y
, mat
, 0.0, max
, rlo
, rhi
, &nlevels
);
457 static void overlap(const char *outfile
, int natoms
,
459 int nvec2
, int *eignr2
, rvec
**eigvec2
,
460 int noutvec
, int *outvec
,
461 const output_env_t oenv
)
467 fprintf(stderr
, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
468 for (i
= 0; i
< noutvec
; i
++)
470 fprintf(stderr
, "%d ", outvec
[i
]+1);
472 fprintf(stderr
, "\n");
474 out
= xvgropen(outfile
, "Subspace overlap",
475 "Eigenvectors of trajectory 2", "Overlap", oenv
);
476 fprintf(out
, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
479 for (x
= 0; x
< nvec2
; x
++)
481 for (v
= 0; v
< noutvec
; v
++)
485 for (i
= 0; i
< natoms
; i
++)
487 inp
+= iprod(eigvec1
[vec
][i
], eigvec2
[x
][i
]);
491 fprintf(out
, "%5d %5.3f\n", eignr2
[x
]+1, overlap
/noutvec
);
497 static void project(const char *trajfile
, t_topology
*top
, int ePBC
, matrix topbox
,
498 const char *projfile
, const char *twodplotfile
,
499 const char *threedplotfile
, const char *filterfile
, int skip
,
500 const char *extremefile
, gmx_bool bExtrAll
, real extreme
,
501 int nextr
, t_atoms
*atoms
, int natoms
, atom_id
*index
,
502 gmx_bool bFit
, rvec
*xref
, int nfit
, atom_id
*ifit
, real
*w_rls
,
503 real
*sqrtm
, rvec
*xav
,
504 int *eignr
, rvec
**eigvec
,
505 int noutvec
, int *outvec
, gmx_bool bSplit
,
506 const output_env_t oenv
)
508 FILE *xvgrout
= NULL
;
509 int nat
, i
, j
, d
, v
, vec
, nfr
, nframes
= 0, snew_size
, frame
;
510 t_trxstatus
*out
= NULL
;
512 int noutvec_extr
, imin
, imax
;
517 real t
, inp
, **inprod
= NULL
, min
= 0, max
= 0;
518 char str
[STRLEN
], str2
[STRLEN
], **ylabel
, *c
;
520 gmx_rmpbc_t gpbc
= NULL
;
526 noutvec_extr
= noutvec
;
536 snew(inprod
, noutvec
+1);
540 fprintf(stderr
, "Writing a filtered trajectory to %s using eigenvectors\n",
542 for (i
= 0; i
< noutvec
; i
++)
544 fprintf(stderr
, "%d ", outvec
[i
]+1);
546 fprintf(stderr
, "\n");
547 out
= open_trx(filterfile
, "w");
552 nat
= read_first_x(oenv
, &status
, trajfile
, &t
, &xread
, box
);
555 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
);
561 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, nat
);
564 for (i
= 0; i
< nat
; i
++)
574 gmx_rmpbc(gpbc
, nat
, box
, xread
);
576 if (nframes
>= snew_size
)
579 for (i
= 0; i
< noutvec
+1; i
++)
581 srenew(inprod
[i
], snew_size
);
584 inprod
[noutvec
][nframes
] = t
;
585 /* calculate x: a fitted struture of the selected atoms */
588 reset_x(nfit
, ifit
, nat
, NULL
, xread
, w_rls
);
589 do_fit(nat
, w_rls
, xref
, xread
);
591 for (i
= 0; i
< natoms
; i
++)
593 copy_rvec(xread
[index
[i
]], x
[i
]);
596 for (v
= 0; v
< noutvec
; v
++)
599 /* calculate (mass-weighted) projection */
601 for (i
= 0; i
< natoms
; i
++)
603 inp
+= (eigvec
[vec
][i
][0]*(x
[i
][0]-xav
[i
][0])+
604 eigvec
[vec
][i
][1]*(x
[i
][1]-xav
[i
][1])+
605 eigvec
[vec
][i
][2]*(x
[i
][2]-xav
[i
][2]))*sqrtm
[i
];
607 inprod
[v
][nframes
] = inp
;
611 for (i
= 0; i
< natoms
; i
++)
613 for (d
= 0; d
< DIM
; d
++)
615 /* misuse xread for output */
616 xread
[index
[i
]][d
] = xav
[i
][d
];
617 for (v
= 0; v
< noutvec
; v
++)
619 xread
[index
[i
]][d
] +=
620 inprod
[v
][nframes
]*eigvec
[outvec
[v
]][i
][d
]/sqrtm
[i
];
624 write_trx(out
, natoms
, index
, atoms
, 0, t
, box
, xread
, NULL
, NULL
);
630 while (read_next_x(oenv
, status
, &t
, xread
, box
));
640 snew(xread
, atoms
->nr
);
645 gmx_rmpbc_done(gpbc
);
651 snew(ylabel
, noutvec
);
652 for (v
= 0; v
< noutvec
; v
++)
654 sprintf(str
, "vec %d", eignr
[outvec
[v
]]+1);
655 ylabel
[v
] = strdup(str
);
657 sprintf(str
, "projection on eigenvectors (%s)", proj_unit
);
658 write_xvgr_graphs(projfile
, noutvec
, 1, str
, NULL
, output_env_get_xvgr_tlabel(oenv
),
659 (const char **)ylabel
,
660 nframes
, inprod
[noutvec
], inprod
, NULL
,
661 output_env_get_time_factor(oenv
), FALSE
, bSplit
, oenv
);
666 sprintf(str
, "projection on eigenvector %d (%s)",
667 eignr
[outvec
[0]]+1, proj_unit
);
668 sprintf(str2
, "projection on eigenvector %d (%s)",
669 eignr
[outvec
[noutvec
-1]]+1, proj_unit
);
670 xvgrout
= xvgropen(twodplotfile
, "2D projection of trajectory", str
, str2
,
672 for (i
= 0; i
< nframes
; i
++)
674 if (bSplit
&& i
> 0 && fabs(inprod
[noutvec
][i
]) < 1e-5)
676 fprintf(xvgrout
, "&\n");
678 fprintf(xvgrout
, "%10.5f %10.5f\n", inprod
[0][i
], inprod
[noutvec
-1][i
]);
680 gmx_ffclose(xvgrout
);
689 char *resnm
, *atnm
, pdbform
[STRLEN
];
695 gmx_fatal(FARGS
, "You have selected less than 3 eigenvectors");
699 bPDB
= fn2ftp(threedplotfile
) == efPDB
;
701 box
[XX
][XX
] = box
[YY
][YY
] = box
[ZZ
][ZZ
] = 1;
703 b4D
= bPDB
&& (noutvec
>= 4);
706 fprintf(stderr
, "You have selected four or more eigenvectors:\n"
707 "fourth eigenvector will be plotted "
708 "in bfactor field of pdb file\n");
709 sprintf(str
, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
710 eignr
[outvec
[0]]+1, eignr
[outvec
[1]]+1,
711 eignr
[outvec
[2]]+1, eignr
[outvec
[3]]+1);
715 sprintf(str
, "3D proj. of traj. on eigenv. %d, %d and %d",
716 eignr
[outvec
[0]]+1, eignr
[outvec
[1]]+1, eignr
[outvec
[2]]+1);
718 init_t_atoms(&atoms
, nframes
, FALSE
);
722 resnm
= strdup("PRJ");
726 fact
= 10000.0/nframes
;
733 for (i
= 0; i
< nframes
; i
++)
735 atoms
.atomname
[i
] = &atnm
;
736 atoms
.atom
[i
].resind
= i
;
737 atoms
.resinfo
[i
].name
= &resnm
;
738 atoms
.resinfo
[i
].nr
= ceil(i
*fact
);
739 atoms
.resinfo
[i
].ic
= ' ';
740 x
[i
][XX
] = inprod
[0][i
];
741 x
[i
][YY
] = inprod
[1][i
];
742 x
[i
][ZZ
] = inprod
[2][i
];
748 if ( ( b4D
|| bSplit
) && bPDB
)
750 strcpy(pdbform
, get_pdbformat());
751 strcat(pdbform
, "%8.4f%8.4f\n");
753 out
= gmx_ffopen(threedplotfile
, "w");
754 fprintf(out
, "HEADER %s\n", str
);
757 fprintf(out
, "REMARK %s\n", "fourth dimension plotted as B-factor");
760 for (i
= 0; i
< atoms
.nr
; i
++)
762 if (j
> 0 && bSplit
&& fabs(inprod
[noutvec
][i
]) < 1e-5)
764 fprintf(out
, "TER\n");
767 fprintf(out
, pdbform
, "ATOM", i
+1, "C", "PRJ", ' ', j
+1,
768 10*x
[i
][XX
], 10*x
[i
][YY
], 10*x
[i
][ZZ
], 1.0, 10*b
[i
]);
771 fprintf(out
, "CONECT%5d%5d\n", i
, i
+1);
775 fprintf(out
, "TER\n");
780 write_sto_conf(threedplotfile
, str
, &atoms
, x
, NULL
, ePBC
, box
);
782 free_t_atoms(&atoms
, FALSE
);
787 snew(pmin
, noutvec_extr
);
788 snew(pmax
, noutvec_extr
);
791 fprintf(stderr
, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
793 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
796 for (v
= 0; v
< noutvec_extr
; v
++)
798 for (i
= 0; i
< nframes
; i
++)
800 if (inprod
[v
][i
] < inprod
[v
][imin
])
804 if (inprod
[v
][i
] > inprod
[v
][imax
])
809 pmin
[v
] = inprod
[v
][imin
];
810 pmax
[v
] = inprod
[v
][imax
];
811 fprintf(stderr
, "%7d %10.6f %10d %10.6f %10d\n",
813 pmin
[v
], imin
, pmax
[v
], imax
);
821 /* build format string for filename: */
822 strcpy(str
, extremefile
); /* copy filename */
823 c
= strrchr(str
, '.'); /* find where extention begins */
824 strcpy(str2
, c
); /* get extention */
825 sprintf(c
, "%%d%s", str2
); /* append '%s' and extention to filename */
826 for (v
= 0; v
< noutvec_extr
; v
++)
828 /* make filename using format string */
829 if (noutvec_extr
== 1)
831 strcpy(str2
, extremefile
);
835 sprintf(str2
, str
, eignr
[outvec
[v
]]+1);
837 fprintf(stderr
, "Writing %d frames along eigenvector %d to %s\n",
838 nextr
, outvec
[v
]+1, str2
);
839 out
= open_trx(str2
, "w");
840 for (frame
= 0; frame
< nextr
; frame
++)
842 if ((extreme
== 0) && (nextr
<= 3))
844 for (i
= 0; i
< natoms
; i
++)
846 atoms
->resinfo
[atoms
->atom
[index
[i
]].resind
].chainid
= 'A' + frame
;
849 for (i
= 0; i
< natoms
; i
++)
851 for (d
= 0; d
< DIM
; d
++)
854 (xav
[i
][d
] + (pmin
[v
]*(nextr
-frame
-1)+pmax
[v
]*frame
)/(nextr
-1)
855 *eigvec
[outvec
[v
]][i
][d
]/sqrtm
[i
]);
858 write_trx(out
, natoms
, index
, atoms
, 0, frame
, topbox
, xread
, NULL
, NULL
);
865 fprintf(stderr
, "\n");
868 static void components(const char *outfile
, int natoms
,
869 int *eignr
, rvec
**eigvec
,
870 int noutvec
, int *outvec
,
871 const output_env_t oenv
)
875 char str
[STRLEN
], **ylabel
;
877 fprintf(stderr
, "Writing eigenvector components to %s\n", outfile
);
879 snew(ylabel
, noutvec
);
882 for (i
= 0; i
< natoms
; i
++)
886 for (g
= 0; g
< noutvec
; g
++)
889 sprintf(str
, "vec %d", eignr
[v
]+1);
890 ylabel
[g
] = strdup(str
);
892 for (s
= 0; s
< 4; s
++)
894 snew(y
[g
][s
], natoms
);
896 for (i
= 0; i
< natoms
; i
++)
898 y
[g
][0][i
] = norm(eigvec
[v
][i
]);
899 for (s
= 0; s
< 3; s
++)
901 y
[g
][s
+1][i
] = eigvec
[v
][i
][s
];
905 write_xvgr_graphs(outfile
, noutvec
, 4, "Eigenvector components",
906 "black: total, red: x, green: y, blue: z",
907 "Atom number", (const char **)ylabel
,
908 natoms
, x
, NULL
, y
, 1, FALSE
, FALSE
, oenv
);
909 fprintf(stderr
, "\n");
912 static void rmsf(const char *outfile
, int natoms
, real
*sqrtm
,
913 int *eignr
, rvec
**eigvec
,
914 int noutvec
, int *outvec
,
915 real
*eigval
, int neig
,
916 const output_env_t oenv
)
920 char str
[STRLEN
], **ylabel
;
922 for (i
= 0; i
< neig
; i
++)
930 fprintf(stderr
, "Writing rmsf to %s\n", outfile
);
932 snew(ylabel
, noutvec
);
935 for (i
= 0; i
< natoms
; i
++)
939 for (g
= 0; g
< noutvec
; g
++)
942 if (eignr
[v
] >= neig
)
944 gmx_fatal(FARGS
, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr
[v
]+1, neig
);
946 sprintf(str
, "vec %d", eignr
[v
]+1);
947 ylabel
[g
] = strdup(str
);
949 for (i
= 0; i
< natoms
; i
++)
951 y
[g
][i
] = sqrt(eigval
[eignr
[v
]]*norm2(eigvec
[v
][i
]))/sqrtm
[i
];
954 write_xvgr_graphs(outfile
, noutvec
, 1, "RMS fluctuation (nm) ", NULL
,
955 "Atom number", (const char **)ylabel
,
956 natoms
, x
, y
, NULL
, 1, TRUE
, FALSE
, oenv
);
957 fprintf(stderr
, "\n");
960 int gmx_anaeig(int argc
, char *argv
[])
962 static const char *desc
[] = {
963 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
964 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
965 "([gmx-nmeig]).[PAR]",
967 "When a trajectory is projected on eigenvectors, all structures are",
968 "fitted to the structure in the eigenvector file, if present, otherwise",
969 "to the structure in the structure file. When no run input file is",
970 "supplied, periodicity will not be taken into account. Most analyses",
971 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
972 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
974 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
975 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
977 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
978 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
980 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
981 "[TT]-first[tt] to [TT]-last[tt].",
982 "The projections of a trajectory on the eigenvectors of its",
983 "covariance matrix are called principal components (pc's).",
984 "It is often useful to check the cosine content of the pc's,",
985 "since the pc's of random diffusion are cosines with the number",
986 "of periods equal to half the pc index.",
987 "The cosine content of the pc's can be calculated with the program",
988 "[gmx-analyze].[PAR]",
990 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
991 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
993 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
994 "three selected eigenvectors.[PAR]",
996 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
997 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
999 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
1000 "on the average structure and interpolate [TT]-nframes[tt] frames",
1001 "between them, or set your own extremes with [TT]-max[tt]. The",
1002 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1003 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1004 "will be written to separate files. Chain identifiers will be added",
1005 "when writing a [TT].pdb[tt] file with two or three structures (you",
1006 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1008 " Overlap calculations between covariance analysis:[BR]",
1009 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1011 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1012 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1013 "in file [TT]-v[tt].[PAR]",
1015 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1016 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1017 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1018 "have been set explicitly.[PAR]",
1020 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1021 "a single number for the overlap between the covariance matrices is",
1022 "generated. The formulas are:[BR]",
1023 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1024 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1025 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1026 "where M1 and M2 are the two covariance matrices and tr is the trace",
1027 "of a matrix. The numbers are proportional to the overlap of the square",
1028 "root of the fluctuations. The normalized overlap is the most useful",
1029 "number, it is 1 for identical matrices and 0 when the sampled",
1030 "subspaces are orthogonal.[PAR]",
1031 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1032 "computed based on the Quasiharmonic approach and based on",
1033 "Schlitter's formula."
1035 static int first
= 1, last
= -1, skip
= 1, nextr
= 2, nskip
= 6;
1036 static real max
= 0.0, temp
= 298.15;
1037 static gmx_bool bSplit
= FALSE
, bEntropy
= FALSE
;
1039 { "-first", FALSE
, etINT
, {&first
},
1040 "First eigenvector for analysis (-1 is select)" },
1041 { "-last", FALSE
, etINT
, {&last
},
1042 "Last eigenvector for analysis (-1 is till the last)" },
1043 { "-skip", FALSE
, etINT
, {&skip
},
1044 "Only analyse every nr-th frame" },
1045 { "-max", FALSE
, etREAL
, {&max
},
1046 "Maximum for projection of the eigenvector on the average structure, "
1047 "max=0 gives the extremes" },
1048 { "-nframes", FALSE
, etINT
, {&nextr
},
1049 "Number of frames for the extremes output" },
1050 { "-split", FALSE
, etBOOL
, {&bSplit
},
1051 "Split eigenvector projections where time is zero" },
1052 { "-entropy", FALSE
, etBOOL
, {&bEntropy
},
1053 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1054 { "-temp", FALSE
, etREAL
, {&temp
},
1055 "Temperature for entropy calculations" },
1056 { "-nevskip", FALSE
, etINT
, {&nskip
},
1057 "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." }
1059 #define NPA asize(pa)
1065 t_atoms
*atoms
= NULL
;
1066 rvec
*xtop
, *xref1
, *xref2
, *xrefp
= NULL
;
1067 gmx_bool bDMR1
, bDMA1
, bDMR2
, bDMA2
;
1068 int nvec1
, nvec2
, *eignr1
= NULL
, *eignr2
= NULL
;
1069 rvec
*x
, *xread
, *xav1
, *xav2
, **eigvec1
= NULL
, **eigvec2
= NULL
;
1071 real xid
, totmass
, *sqrtm
, *w_rls
, t
, lambda
;
1074 const char *indexfile
;
1077 int nout
, *iout
, noutvec
, *outvec
, nfit
;
1078 atom_id
*index
, *ifit
;
1079 const char *VecFile
, *Vec2File
, *topfile
;
1080 const char *EigFile
, *Eig2File
;
1081 const char *CompFile
, *RmsfFile
, *ProjOnVecFile
;
1082 const char *TwoDPlotFile
, *ThreeDPlotFile
;
1083 const char *FilterFile
, *ExtremeFile
;
1084 const char *OverlapFile
, *InpMatFile
;
1085 gmx_bool bFit1
, bFit2
, bM
, bIndex
, bTPS
, bTop
, bVec2
, bProj
;
1086 gmx_bool bFirstToLast
, bFirstLastSet
, bTraj
, bCompare
, bPDB3D
;
1087 real
*eigval1
= NULL
, *eigval2
= NULL
;
1094 { efTRN
, "-v", "eigenvec", ffREAD
},
1095 { efTRN
, "-v2", "eigenvec2", ffOPTRD
},
1096 { efTRX
, "-f", NULL
, ffOPTRD
},
1097 { efTPS
, NULL
, NULL
, ffOPTRD
},
1098 { efNDX
, NULL
, NULL
, ffOPTRD
},
1099 { efXVG
, "-eig", "eigenval", ffOPTRD
},
1100 { efXVG
, "-eig2", "eigenval2", ffOPTRD
},
1101 { efXVG
, "-comp", "eigcomp", ffOPTWR
},
1102 { efXVG
, "-rmsf", "eigrmsf", ffOPTWR
},
1103 { efXVG
, "-proj", "proj", ffOPTWR
},
1104 { efXVG
, "-2d", "2dproj", ffOPTWR
},
1105 { efSTO
, "-3d", "3dproj.pdb", ffOPTWR
},
1106 { efTRX
, "-filt", "filtered", ffOPTWR
},
1107 { efTRX
, "-extr", "extreme.pdb", ffOPTWR
},
1108 { efXVG
, "-over", "overlap", ffOPTWR
},
1109 { efXPM
, "-inpr", "inprod", ffOPTWR
}
1111 #define NFILE asize(fnm)
1113 if (!parse_common_args(&argc
, argv
,
1114 PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_CAN_VIEW
| PCA_BE_NICE
,
1115 NFILE
, fnm
, NPA
, pa
, asize(desc
), desc
, 0, NULL
, &oenv
))
1120 indexfile
= ftp2fn_null(efNDX
, NFILE
, fnm
);
1122 VecFile
= opt2fn("-v", NFILE
, fnm
);
1123 Vec2File
= opt2fn_null("-v2", NFILE
, fnm
);
1124 topfile
= ftp2fn(efTPS
, NFILE
, fnm
);
1125 EigFile
= opt2fn_null("-eig", NFILE
, fnm
);
1126 Eig2File
= opt2fn_null("-eig2", NFILE
, fnm
);
1127 CompFile
= opt2fn_null("-comp", NFILE
, fnm
);
1128 RmsfFile
= opt2fn_null("-rmsf", NFILE
, fnm
);
1129 ProjOnVecFile
= opt2fn_null("-proj", NFILE
, fnm
);
1130 TwoDPlotFile
= opt2fn_null("-2d", NFILE
, fnm
);
1131 ThreeDPlotFile
= opt2fn_null("-3d", NFILE
, fnm
);
1132 FilterFile
= opt2fn_null("-filt", NFILE
, fnm
);
1133 ExtremeFile
= opt2fn_null("-extr", NFILE
, fnm
);
1134 OverlapFile
= opt2fn_null("-over", NFILE
, fnm
);
1135 InpMatFile
= ftp2fn_null(efXPM
, NFILE
, fnm
);
1137 bTop
= fn2bTPX(topfile
);
1138 bProj
= ProjOnVecFile
|| TwoDPlotFile
|| ThreeDPlotFile
1139 || FilterFile
|| ExtremeFile
;
1141 opt2parg_bSet("-first", NPA
, pa
) && opt2parg_bSet("-last", NPA
, pa
);
1142 bFirstToLast
= CompFile
|| RmsfFile
|| ProjOnVecFile
|| FilterFile
||
1143 OverlapFile
|| ((ExtremeFile
|| InpMatFile
) && bFirstLastSet
);
1144 bVec2
= Vec2File
|| OverlapFile
|| InpMatFile
;
1145 bM
= RmsfFile
|| bProj
;
1146 bTraj
= ProjOnVecFile
|| FilterFile
|| (ExtremeFile
&& (max
== 0))
1147 || TwoDPlotFile
|| ThreeDPlotFile
;
1148 bIndex
= bM
|| bProj
;
1149 bTPS
= ftp2bSet(efTPS
, NFILE
, fnm
) || bM
|| bTraj
||
1150 FilterFile
|| (bIndex
&& indexfile
);
1151 bCompare
= Vec2File
|| Eig2File
;
1152 bPDB3D
= fn2ftp(ThreeDPlotFile
) == efPDB
;
1154 read_eigenvectors(VecFile
, &natoms
, &bFit1
,
1155 &xref1
, &bDMR1
, &xav1
, &bDMA1
,
1156 &nvec1
, &eignr1
, &eigvec1
, &eigval1
);
1159 /* Overwrite eigenvalues from separate files if the user provides them */
1160 if (EigFile
!= NULL
)
1162 int neig_tmp
= read_xvg(EigFile
, &xvgdata
, &i
);
1163 if (neig_tmp
!= neig1
)
1165 fprintf(stderr
, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1
, natoms
);
1168 srenew(eigval1
, neig1
);
1169 for (j
= 0; j
< neig1
; j
++)
1171 real tmp
= eigval1
[j
];
1172 eigval1
[j
] = xvgdata
[1][j
];
1173 if (debug
&& (eigval1
[j
] != tmp
))
1175 fprintf(debug
, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1176 j
, tmp
, eigval1
[j
]);
1179 for (j
= 0; j
< i
; j
++)
1184 fprintf(stderr
, "Read %d eigenvalues from %s\n", neig1
, EigFile
);
1191 gmx_fatal(FARGS
, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1193 calc_entropy_qh(stdout
, neig1
, eigval1
, temp
, nskip
);
1194 calc_entropy_schlitter(stdout
, neig1
, nskip
, eigval1
, temp
);
1201 gmx_fatal(FARGS
, "Need a second eigenvector file to do this analysis.");
1203 read_eigenvectors(Vec2File
, &neig2
, &bFit2
,
1204 &xref2
, &bDMR2
, &xav2
, &bDMA2
, &nvec2
, &eignr2
, &eigvec2
, &eigval2
);
1209 gmx_fatal(FARGS
, "Dimensions in the eigenvector files don't match");
1213 if (Eig2File
!= NULL
)
1215 neig2
= read_xvg(Eig2File
, &xvgdata
, &i
);
1216 srenew(eigval2
, neig2
);
1217 for (j
= 0; j
< neig2
; j
++)
1219 eigval2
[j
] = xvgdata
[1][j
];
1221 for (j
= 0; j
< i
; j
++)
1226 fprintf(stderr
, "Read %d eigenvalues from %s\n", neig2
, Eig2File
);
1230 if ((!bFit1
|| xref1
) && !bDMR1
&& !bDMA1
)
1234 if ((xref1
== NULL
) && (bM
|| bTraj
))
1250 bTop
= read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
),
1251 title
, &top
, &ePBC
, &xtop
, NULL
, topbox
, bM
);
1253 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, atoms
->nr
);
1254 gmx_rmpbc(gpbc
, atoms
->nr
, topbox
, xtop
);
1255 /* Fitting is only required for the projection */
1260 printf("\nNote: the structure in %s should be the same\n"
1261 " as the one used for the fit in g_covar\n", topfile
);
1263 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1264 get_index(atoms
, indexfile
, 1, &nfit
, &ifit
, &grpname
);
1266 snew(w_rls
, atoms
->nr
);
1267 for (i
= 0; (i
< nfit
); i
++)
1271 w_rls
[ifit
[i
]] = atoms
->atom
[ifit
[i
]].m
;
1275 w_rls
[ifit
[i
]] = 1.0;
1279 snew(xrefp
, atoms
->nr
);
1282 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1285 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
);
1287 for (i
= 0; (i
< nfit
); i
++)
1289 copy_rvec(xref1
[i
], xrefp
[ifit
[i
]]);
1294 /* The top coordinates are the fitting reference */
1295 for (i
= 0; (i
< nfit
); i
++)
1297 copy_rvec(xtop
[ifit
[i
]], xrefp
[ifit
[i
]]);
1299 reset_x(nfit
, ifit
, atoms
->nr
, NULL
, xrefp
, w_rls
);
1302 gmx_rmpbc_done(gpbc
);
1307 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms
);
1308 get_index(atoms
, indexfile
, 1, &i
, &index
, &grpname
);
1311 gmx_fatal(FARGS
, "you selected a group with %d elements instead of %d", i
, natoms
);
1316 snew(sqrtm
, natoms
);
1319 proj_unit
= "u\\S1/2\\Nnm";
1320 for (i
= 0; (i
< natoms
); i
++)
1322 sqrtm
[i
] = sqrt(atoms
->atom
[index
[i
]].m
);
1328 for (i
= 0; (i
< natoms
); i
++)
1338 for (i
= 0; (i
< natoms
); i
++)
1340 for (d
= 0; (d
< DIM
); d
++)
1342 t
+= sqr((xav1
[i
][d
]-xav2
[i
][d
])*sqrtm
[i
]);
1343 totmass
+= sqr(sqrtm
[i
]);
1346 fprintf(stdout
, "RMSD (without fit) between the two average structures:"
1347 " %.3f (nm)\n\n", sqrt(t
/totmass
));
1358 /* make an index from first to last */
1359 nout
= last
-first
+1;
1361 for (i
= 0; i
< nout
; i
++)
1363 iout
[i
] = first
-1+i
;
1366 else if (ThreeDPlotFile
)
1368 /* make an index of first+(0,1,2) and last */
1369 nout
= bPDB3D
? 4 : 3;
1370 nout
= min(last
-first
+1, nout
);
1378 iout
[nout
-1] = last
-1;
1382 /* make an index of first and last */
1391 printf("Select eigenvectors for output, end your selection with 0\n");
1398 srenew(iout
, nout
+1);
1399 if (1 != scanf("%d", &iout
[nout
]))
1401 gmx_fatal(FARGS
, "Error reading user input");
1405 while (iout
[nout
] >= 0);
1409 /* make an index of the eigenvectors which are present */
1412 for (i
= 0; i
< nout
; i
++)
1415 while ((j
< nvec1
) && (eignr1
[j
] != iout
[i
]))
1419 if ((j
< nvec1
) && (eignr1
[j
] == iout
[i
]))
1421 outvec
[noutvec
] = j
;
1425 fprintf(stderr
, "%d eigenvectors selected for output", noutvec
);
1428 fprintf(stderr
, ":");
1429 for (j
= 0; j
< noutvec
; j
++)
1431 fprintf(stderr
, " %d", eignr1
[outvec
[j
]]+1);
1434 fprintf(stderr
, "\n");
1438 components(CompFile
, natoms
, eignr1
, eigvec1
, noutvec
, outvec
, oenv
);
1443 rmsf(RmsfFile
, natoms
, sqrtm
, eignr1
, eigvec1
, noutvec
, outvec
, eigval1
,
1449 project(bTraj
? opt2fn("-f", NFILE
, fnm
) : NULL
,
1450 bTop
? &top
: NULL
, ePBC
, topbox
,
1451 ProjOnVecFile
, TwoDPlotFile
, ThreeDPlotFile
, FilterFile
, skip
,
1452 ExtremeFile
, bFirstLastSet
, max
, nextr
, atoms
, natoms
, index
,
1453 bFit1
, xrefp
, nfit
, ifit
, w_rls
,
1454 sqrtm
, xav1
, eignr1
, eigvec1
, noutvec
, outvec
, bSplit
,
1460 overlap(OverlapFile
, natoms
,
1461 eigvec1
, nvec2
, eignr2
, eigvec2
, noutvec
, outvec
, oenv
);
1466 inprod_matrix(InpMatFile
, natoms
,
1467 nvec1
, eignr1
, eigvec1
, nvec2
, eignr2
, eigvec2
,
1468 bFirstLastSet
, noutvec
, outvec
);
1473 compare(natoms
, nvec1
, eigvec1
, nvec2
, eigvec2
, eigval1
, neig1
, eigval2
, neig2
);
1477 if (!CompFile
&& !bProj
&& !OverlapFile
&& !InpMatFile
&&
1478 !bCompare
&& !bEntropy
)
1480 fprintf(stderr
, "\nIf you want some output,"
1481 " set one (or two or ...) of the output file options\n");
1485 view_all(oenv
, NFILE
, fnm
);