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.
47 #include "gromacs/commandline/pargs.h"
50 #include "gmx_fatal.h"
51 #include "gromacs/fileio/futil.h"
55 #include "gromacs/fileio/matio.h"
56 #include "gromacs/fileio/tpxio.h"
57 #include "gromacs/fileio/trxio.h"
63 static void norm_princ(t_atoms *atoms, int isize, atom_id *index, int natoms,
69 /* equalize principal components: */
70 /* orient principal axes, get principal components */
71 orient_princ(atoms, isize, index, natoms, x, NULL, princ);
73 /* calc our own principal components */
75 for (m = 0; m < DIM; m++)
77 for (i = 0; i < isize; i++)
79 vec[m] += sqr(x[index[i]][m]);
81 vec[m] = sqrt(vec[m] / isize);
82 /* calculate scaling constants */
83 vec[m] = 1 / (sqrt(3) * vec[m]);
86 /* scale coordinates */
87 for (i = 0; i < natoms; i++)
89 for (m = 0; m < DIM; m++)
96 int gmx_rms(int argc, char *argv[])
100 "[THISMODULE] compares two structures by computing the root mean square",
101 "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
102 "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
103 "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
104 "This is selected by [TT]-what[tt].[PAR]"
106 "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
107 "reference structure. The reference structure",
108 "is taken from the structure file ([TT]-s[tt]).[PAR]",
110 "With option [TT]-mir[tt] also a comparison with the mirror image of",
111 "the reference structure is calculated.",
112 "This is useful as a reference for 'significant' values, see",
113 "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
115 "Option [TT]-prev[tt] produces the comparison with a previous frame",
116 "the specified number of frames ago.[PAR]",
118 "Option [TT]-m[tt] produces a matrix in [TT].xpm[tt] format of",
119 "comparison values of each structure in the trajectory with respect to",
120 "each other structure. This file can be visualized with for instance",
121 "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]",
123 "Option [TT]-fit[tt] controls the least-squares fitting of",
124 "the structures on top of each other: complete fit (rotation and",
125 "translation), translation only, or no fitting at all.[PAR]",
127 "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
128 "If you select the option (default) and ",
129 "supply a valid [TT].tpr[tt] file masses will be taken from there, ",
130 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
131 "[TT]GMXLIB[tt]. This is fine for proteins, but not",
132 "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
133 "is assigned to unknown atoms. You can check whether this happend by",
134 "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
136 "With [TT]-f2[tt], the 'other structures' are taken from a second",
137 "trajectory, this generates a comparison matrix of one trajectory",
138 "versus the other.[PAR]",
140 "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
142 "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
143 "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
144 "comparison group are considered."
146 static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
147 static gmx_bool bDeltaLog = FALSE;
148 static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
149 static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1,
150 bond_user_min = -1, delta_maxy = 0.0;
151 /* strings and things for selecting difference method */
154 ewSel, ewRMSD, ewRho, ewRhoSc, ewNR
157 const char *what[ewNR + 1] =
158 { NULL, "rmsd", "rho", "rhosc", NULL };
159 const char *whatname[ewNR] =
160 { NULL, "RMSD", "Rho", "Rho sc" };
161 const char *whatlabel[ewNR] =
162 { NULL, "RMSD (nm)", "Rho", "Rho sc" };
163 const char *whatxvgname[ewNR] =
164 { NULL, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
165 const char *whatxvglabel[ewNR] =
166 { NULL, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
167 /* strings and things for fitting methods */
170 efSel, efFit, efReset, efNone, efNR
173 const char *fit[efNR + 1] =
174 { NULL, "rot+trans", "translation", "none", NULL };
175 const char *fitgraphlabel[efNR + 1] =
176 { NULL, "lsq fit", "translational fit", "no fit" };
178 static gmx_bool bMassWeighted = TRUE;
181 { "-what", FALSE, etENUM,
182 { what }, "Structural difference measure" },
183 { "-pbc", FALSE, etBOOL,
184 { &bPBC }, "PBC check" },
185 { "-fit", FALSE, etENUM,
186 { fit }, "Fit to reference structure" },
187 { "-prev", FALSE, etINT,
188 { &prev }, "Compare with previous frame" },
189 { "-split", FALSE, etBOOL,
190 { &bSplit }, "Split graph where time is zero" },
191 { "-fitall", FALSE, etBOOL,
192 { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
193 { "-skip", FALSE, etINT,
194 { &freq }, "Only write every nr-th frame to matrix" },
195 { "-skip2", FALSE, etINT,
196 { &freq2 }, "Only write every nr-th frame to matrix" },
197 { "-max", FALSE, etREAL,
198 { &rmsd_user_max }, "Maximum level in comparison matrix" },
199 { "-min", FALSE, etREAL,
200 { &rmsd_user_min }, "Minimum level in comparison matrix" },
201 { "-bmax", FALSE, etREAL,
202 { &bond_user_max }, "Maximum level in bond angle matrix" },
203 { "-bmin", FALSE, etREAL,
204 { &bond_user_min }, "Minimum level in bond angle matrix" },
205 { "-mw", FALSE, etBOOL,
206 { &bMassWeighted }, "Use mass weighting for superposition" },
207 { "-nlevels", FALSE, etINT,
208 { &nlevels }, "Number of levels in the matrices" },
209 { "-ng", FALSE, etINT,
210 { &nrms }, "Number of groups to compute RMS between" },
211 { "-dlog", FALSE, etBOOL,
213 "HIDDENUse a log x-axis in the delta t matrix" },
214 { "-dmax", FALSE, etREAL,
215 { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
216 { "-aver", FALSE, etINT,
218 "HIDDENAverage over this distance in the RMSD matrix" }
220 int natoms_trx, natoms_trx2, natoms;
221 int i, j, k, m, teller, teller2, tel_mat, tel_mat2;
223 int maxframe = NFRAME, maxframe2 = NFRAME;
224 real t, *w_rls, *w_rms, *w_rls_m = NULL, *w_rms_m = NULL;
225 gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
226 gmx_bool bFit, bReset;
229 t_iatom *iatom = NULL;
232 rvec *x, *xp, *xm = NULL, **mat_x = NULL, **mat_x2, *mat_x2_j = NULL, vec1,
235 char buf[256], buf2[256];
238 real rlstot = 0, **rls, **rlsm = NULL, *time, *time2, *rlsnorm = NULL,
239 **rmsd_mat = NULL, **bond_mat = NULL, *axis, *axis2, *del_xaxis,
240 *del_yaxis, rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
241 real **rmsdav_mat = NULL, av_tot, weight, weight_tot;
242 real **delta = NULL, delta_max, delta_scalex = 0, delta_scaley = 0,
244 int delta_xsize = 0, del_lev = 100, mx, my, abs_my;
245 gmx_bool bA1, bA2, bPrev, bTop, *bInMat = NULL;
246 int ifit, *irms, ibond = 0, *ind_bond1 = NULL, *ind_bond2 = NULL, n_ind_m =
248 atom_id *ind_fit, **ind_rms, *ind_m = NULL, *rev_ind_m = NULL, *ind_rms_m =
250 char *gn_fit, **gn_rms;
253 gmx_rmpbc_t gpbc = NULL;
257 { efTPS, NULL, NULL, ffREAD },
258 { efTRX, "-f", NULL, ffREAD },
259 { efTRX, "-f2", NULL, ffOPTRD },
260 { efNDX, NULL, NULL, ffOPTRD },
261 { efXVG, NULL, "rmsd", ffWRITE },
262 { efXVG, "-mir", "rmsdmir", ffOPTWR },
263 { efXVG, "-a", "avgrp", ffOPTWR },
264 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
265 { efXPM, "-m", "rmsd", ffOPTWR },
266 { efDAT, "-bin", "rmsd", ffOPTWR },
267 { efXPM, "-bm", "bond", ffOPTWR }
269 #define NFILE asize(fnm)
271 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW
272 | PCA_BE_NICE, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
277 /* parse enumerated options: */
279 if (ewhat == ewRho || ewhat == ewRhoSc)
281 please_cite(stdout, "Maiorov95");
284 bFit = efit == efFit;
285 bReset = efit == efReset;
288 bReset = TRUE; /* for fit, reset *must* be set */
295 /* mark active cmdline options */
296 bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
297 bFile2 = opt2bSet("-f2", NFILE, fnm);
298 bMat = opt2bSet("-m", NFILE, fnm);
299 bBond = opt2bSet("-bm", NFILE, fnm);
300 bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
301 * your RMSD matrix (hidden option */
302 bNorm = opt2bSet("-a", NFILE, fnm);
303 bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
306 fprintf(stderr, "The number of frames to skip is <= 0. "
307 "Writing out all frames.\n\n");
314 else if (bFile2 && freq2 <= 0)
317 "The number of frames to skip in second trajectory is <= 0.\n"
318 " Writing out all frames.\n\n");
328 fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
332 if (bFile2 && !bMat && !bBond)
336 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
337 " will not read from %s\n", opt2fn("-f2", NFILE,
348 "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
349 " will not read from %s\n", opt2fn("-f2",
355 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xp,
357 snew(w_rls, top.atoms.nr);
358 snew(w_rms, top.atoms.nr);
363 "WARNING: Need a run input file for bond angle matrix,\n"
364 " will not calculate bond angle matrix.\n");
370 fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares"
372 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit,
382 if (bFit && ifit < 3)
384 gmx_fatal(FARGS, "Need >= 3 points to fit!\n" );
388 for (i = 0; i < ifit; i++)
392 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
396 w_rls[ind_fit[i]] = 1;
398 bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
402 fprintf(stderr, "All masses in the fit group are 0, using masses of 1\n");
403 for (i = 0; i < ifit; i++)
405 w_rls[ind_fit[i]] = 1;
419 fprintf(stderr, "Select group%s for %s calculation\n",
420 (nrms > 1) ? "s" : "", whatname[ewhat]);
421 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm),
422 nrms, irms, ind_rms, gn_rms);
426 snew(rlsnorm, irms[0]);
429 for (j = 0; j < nrms; j++)
431 snew(rls[j], maxframe);
436 for (j = 0; j < nrms; j++)
438 snew(rlsm[j], maxframe);
441 snew(time, maxframe);
442 for (j = 0; j < nrms; j++)
445 for (i = 0; i < irms[j]; i++)
449 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
453 w_rms[ind_rms[j][i]] = 1.0;
455 bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
459 fprintf(stderr, "All masses in group %d are 0, using masses of 1\n", j);
460 for (i = 0; i < irms[j]; i++)
462 w_rms[ind_rms[j][i]] = 1;
466 /* Prepare reference frame */
469 gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr);
470 gmx_rmpbc(gpbc, top.atoms.nr, box, xp);
474 reset_x(ifit, ind_fit, top.atoms.nr, NULL, xp, w_rls);
478 /* generate reference structure mirror image: */
479 snew(xm, top.atoms.nr);
480 for (i = 0; i < top.atoms.nr; i++)
482 copy_rvec(xp[i], xm[i]);
483 xm[i][XX] = -xm[i][XX];
486 if (ewhat == ewRhoSc)
488 norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
491 /* read first frame */
492 natoms_trx = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
493 if (natoms_trx != top.atoms.nr)
496 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
497 top.atoms.nr, natoms_trx);
499 natoms = min(top.atoms.nr, natoms_trx);
500 if (bMat || bBond || bPrev)
506 /* With -prev we use all atoms for simplicity */
511 /* Check which atoms we need (fit/rms) */
512 snew(bInMat, natoms);
513 for (i = 0; i < ifit; i++)
515 bInMat[ind_fit[i]] = TRUE;
518 for (i = 0; i < irms[0]; i++)
520 if (!bInMat[ind_rms[0][i]])
522 bInMat[ind_rms[0][i]] = TRUE;
527 /* Make an index of needed atoms */
528 snew(ind_m, n_ind_m);
529 snew(rev_ind_m, natoms);
531 for (i = 0; i < natoms; i++)
533 if (bPrev || bInMat[i])
540 snew(w_rls_m, n_ind_m);
541 snew(ind_rms_m, irms[0]);
542 snew(w_rms_m, n_ind_m);
543 for (i = 0; i < ifit; i++)
545 w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
547 for (i = 0; i < irms[0]; i++)
549 ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
550 w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
558 for (k = 0; k < F_NRE; k++)
562 iatom = top.idef.il[k].iatoms;
563 ncons += top.idef.il[k].nr/3;
566 fprintf(stderr, "Found %d bonds in topology\n", ncons);
567 snew(ind_bond1, ncons);
568 snew(ind_bond2, ncons);
570 for (k = 0; k < F_NRE; k++)
574 iatom = top.idef.il[k].iatoms;
575 ncons = top.idef.il[k].nr/3;
576 for (i = 0; i < ncons; i++)
580 for (j = 0; j < irms[0]; j++)
582 if (iatom[3*i+1] == ind_rms[0][j])
586 if (iatom[3*i+2] == ind_rms[0][j])
593 ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
594 ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
600 fprintf(stderr, "Using %d bonds for bond angle matrix\n", ibond);
603 gmx_fatal(FARGS, "0 bonds found");
607 /* start looping over frames: */
614 gmx_rmpbc(gpbc, natoms, box, x);
619 reset_x(ifit, ind_fit, natoms, NULL, x, w_rls);
621 if (ewhat == ewRhoSc)
623 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
628 /*do the least squares fit to original structure*/
629 do_fit(natoms, w_rls, xp, x);
632 if (teller % freq == 0)
634 /* keep frame for matrix calculation */
635 if (bMat || bBond || bPrev)
637 if (tel_mat >= NFRAME)
639 srenew(mat_x, tel_mat+1);
641 snew(mat_x[tel_mat], n_ind_m);
642 for (i = 0; i < n_ind_m; i++)
644 copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]);
650 /*calculate energy of root_least_squares*/
658 for (i = 0; i < n_ind_m; i++)
660 copy_rvec(mat_x[j][i], xp[ind_m[i]]);
664 reset_x(ifit, ind_fit, natoms, NULL, xp, w_rls);
668 do_fit(natoms, w_rls, x, xp);
671 for (j = 0; (j < nrms); j++)
674 calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
678 for (j = 0; (j < irms[0]); j++)
681 calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
689 /*do the least squares fit to mirror of original structure*/
690 do_fit(natoms, w_rls, xm, x);
693 for (j = 0; j < nrms; j++)
696 calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
699 time[teller] = output_env_conv_time(oenv, t);
702 if (teller >= maxframe)
705 srenew(time, maxframe);
706 for (j = 0; (j < nrms); j++)
708 srenew(rls[j], maxframe);
712 for (j = 0; (j < nrms); j++)
714 srenew(rlsm[j], maxframe);
719 while (read_next_x(oenv, status, &t, x, box));
724 snew(time2, maxframe2);
726 fprintf(stderr, "\nWill read second trajectory file\n");
727 snew(mat_x2, NFRAME);
729 read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
730 if (natoms_trx2 != natoms_trx)
733 "Second trajectory (%d atoms) does not match the first one"
734 " (%d atoms)", natoms_trx2, natoms_trx);
742 gmx_rmpbc(gpbc, natoms, box, x);
747 reset_x(ifit, ind_fit, natoms, NULL, x, w_rls);
749 if (ewhat == ewRhoSc)
751 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
756 /*do the least squares fit to original structure*/
757 do_fit(natoms, w_rls, xp, x);
760 if (teller2 % freq2 == 0)
762 /* keep frame for matrix calculation */
765 if (tel_mat2 >= NFRAME)
767 srenew(mat_x2, tel_mat2+1);
769 snew(mat_x2[tel_mat2], n_ind_m);
770 for (i = 0; i < n_ind_m; i++)
772 copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]);
778 time2[teller2] = output_env_conv_time(oenv, t);
781 if (teller2 >= maxframe2)
784 srenew(time2, maxframe2);
787 while (read_next_x(oenv, status, &t, x, box));
797 gmx_rmpbc_done(gpbc);
801 /* calculate RMS matrix */
802 fprintf(stderr, "\n");
805 fprintf(stderr, "Building %s matrix, %dx%d elements\n",
806 whatname[ewhat], tel_mat, tel_mat2);
807 snew(rmsd_mat, tel_mat);
811 fprintf(stderr, "Building bond angle matrix, %dx%d elements\n",
813 snew(bond_mat, tel_mat);
816 snew(axis2, tel_mat2);
829 for (j = 0; j < tel_mat2; j++)
831 axis2[j] = time2[freq2*j];
837 delta_scalex = 8.0/log(2.0);
838 delta_xsize = (int)(log(tel_mat/2)*delta_scalex+0.5)+1;
842 delta_xsize = tel_mat/2;
844 delta_scaley = 1.0/delta_maxy;
845 snew(delta, delta_xsize);
846 for (j = 0; j < delta_xsize; j++)
848 snew(delta[j], del_lev+1);
852 snew(rmsdav_mat, tel_mat);
853 for (j = 0; j < tel_mat; j++)
855 snew(rmsdav_mat[j], tel_mat);
862 snew(mat_x2_j, natoms);
864 for (i = 0; i < tel_mat; i++)
866 axis[i] = time[freq*i];
867 fprintf(stderr, "\r element %5d; time %5.2f ", i, axis[i]);
870 snew(rmsd_mat[i], tel_mat2);
874 snew(bond_mat[i], tel_mat2);
876 for (j = 0; j < tel_mat2; j++)
880 for (k = 0; k < n_ind_m; k++)
882 copy_rvec(mat_x2[j][k], mat_x2_j[k]);
884 do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j);
888 mat_x2_j = mat_x2[j];
892 if (bFile2 || (i < j))
895 calc_similar_ind(ewhat != ewRMSD, irms[0], ind_rms_m,
896 w_rms_m, mat_x[i], mat_x2_j);
897 if (rmsd_mat[i][j] > rmsd_max)
899 rmsd_max = rmsd_mat[i][j];
901 if (rmsd_mat[i][j] < rmsd_min)
903 rmsd_min = rmsd_mat[i][j];
905 rmsd_avg += rmsd_mat[i][j];
909 rmsd_mat[i][j] = rmsd_mat[j][i];
914 if (bFile2 || (i <= j))
917 for (m = 0; m < ibond; m++)
919 rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1);
920 rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2);
921 ang += acos(cos_angle(vec1, vec2));
923 bond_mat[i][j] = ang*180.0/(M_PI*ibond);
924 if (bond_mat[i][j] > bond_max)
926 bond_max = bond_mat[i][j];
928 if (bond_mat[i][j] < bond_min)
930 bond_min = bond_mat[i][j];
935 bond_mat[i][j] = bond_mat[j][i];
942 rmsd_avg /= tel_mat*tel_mat2;
946 rmsd_avg /= tel_mat*(tel_mat - 1)/2;
948 if (bMat && (avl > 0))
953 for (j = 0; j < tel_mat-1; j++)
955 for (i = j+1; i < tel_mat; i++)
959 for (my = -avl; my <= avl; my++)
961 if ((j+my >= 0) && (j+my < tel_mat))
964 for (mx = -avl; mx <= avl; mx++)
966 if ((i+mx >= 0) && (i+mx < tel_mat))
968 weight = (real)(avl+1-max(abs(mx), abs_my));
969 av_tot += weight*rmsd_mat[i+mx][j+my];
970 weight_tot += weight;
975 rmsdav_mat[i][j] = av_tot/weight_tot;
976 rmsdav_mat[j][i] = rmsdav_mat[i][j];
977 if (rmsdav_mat[i][j] > rmsd_max)
979 rmsd_max = rmsdav_mat[i][j];
983 rmsd_mat = rmsdav_mat;
988 fprintf(stderr, "\n%s: Min %f, Max %f, Avg %f\n",
989 whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg);
990 rlo.r = 1; rlo.g = 1; rlo.b = 1;
991 rhi.r = 0; rhi.g = 0; rhi.b = 0;
992 if (rmsd_user_max != -1)
994 rmsd_max = rmsd_user_max;
996 if (rmsd_user_min != -1)
998 rmsd_min = rmsd_user_min;
1000 if ((rmsd_user_max != -1) || (rmsd_user_min != -1))
1002 fprintf(stderr, "Min and Max value set to resp. %f and %f\n",
1003 rmsd_min, rmsd_max);
1005 sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]);
1006 write_xpm(opt2FILE("-m", NFILE, fnm, "w"), 0, buf, whatlabel[ewhat],
1007 output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
1008 axis, axis2, rmsd_mat, rmsd_min, rmsd_max, rlo, rhi, &nlevels);
1009 /* Print the distribution of RMSD values */
1010 if (opt2bSet("-dist", NFILE, fnm))
1012 low_rmsd_dist(opt2fn("-dist", NFILE, fnm), rmsd_max, tel_mat, rmsd_mat, oenv);
1017 snew(delta_tot, delta_xsize);
1018 for (j = 0; j < tel_mat-1; j++)
1020 for (i = j+1; i < tel_mat; i++)
1027 mx = (int)(log(mx)*delta_scalex+0.5);
1029 my = (int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
1030 delta_tot[mx] += 1.0;
1031 if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy))
1033 delta[mx][my] += 1.0;
1039 for (i = 0; i < delta_xsize; i++)
1041 if (delta_tot[i] > 0.0)
1043 delta_tot[i] = 1.0/delta_tot[i];
1044 for (j = 0; j <= del_lev; j++)
1046 delta[i][j] *= delta_tot[i];
1047 if (delta[i][j] > delta_max)
1049 delta_max = delta[i][j];
1054 fprintf(stderr, "Maximum in delta matrix: %f\n", delta_max);
1055 snew(del_xaxis, delta_xsize);
1056 snew(del_yaxis, del_lev+1);
1057 for (i = 0; i < delta_xsize; i++)
1059 del_xaxis[i] = axis[i]-axis[0];
1061 for (i = 0; i < del_lev+1; i++)
1063 del_yaxis[i] = delta_maxy*i/del_lev;
1065 sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]);
1066 fp = gmx_ffopen("delta.xpm", "w");
1067 write_xpm(fp, 0, buf, "density", output_env_get_time_label(oenv), whatlabel[ewhat],
1068 delta_xsize, del_lev+1, del_xaxis, del_yaxis,
1069 delta, 0.0, delta_max, rlo, rhi, &nlevels);
1072 if (opt2bSet("-bin", NFILE, fnm))
1074 /* NB: File must be binary if we use fwrite */
1075 fp = ftp2FILE(efDAT, NFILE, fnm, "wb");
1076 for (i = 0; i < tel_mat; i++)
1078 if (fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp) != tel_mat2)
1080 gmx_fatal(FARGS, "Error writing to output file");
1088 fprintf(stderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1089 if (bond_user_max != -1)
1091 bond_max = bond_user_max;
1093 if (bond_user_min != -1)
1095 bond_min = bond_user_min;
1097 if ((bond_user_max != -1) || (bond_user_min != -1))
1099 fprintf(stderr, "Bond angle Min and Max set to:\n"
1100 "Min. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1102 rlo.r = 1; rlo.g = 1; rlo.b = 1;
1103 rhi.r = 0; rhi.g = 0; rhi.b = 0;
1104 sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);
1105 write_xpm(opt2FILE("-bm", NFILE, fnm, "w"), 0, buf, "degrees",
1106 output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2,
1107 axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels);
1111 bAv = opt2bSet("-a", NFILE, fnm);
1113 /* Write the RMSD's to file */
1116 sprintf(buf, "%s", whatxvgname[ewhat]);
1120 sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat],
1121 time[prev*freq]-time[0], output_env_get_time_label(oenv));
1123 fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
1124 whatxvglabel[ewhat], oenv);
1125 if (output_env_get_print_xvgr_codes(oenv))
1127 fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n",
1128 (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit],
1129 bFit ? " to " : "", bFit ? gn_fit : "");
1133 xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
1135 for (i = 0; (i < teller); i++)
1137 if (bSplit && i > 0 &&
1138 abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5)
1142 fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]);
1143 for (j = 0; (j < nrms); j++)
1145 fprintf(fp, " %12.7f", rls[j][i]);
1148 rlstot += rls[j][i];
1157 /* Write the mirror RMSD's to file */
1158 sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
1159 sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
1160 fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
1164 if (output_env_get_print_xvgr_codes(oenv))
1166 fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
1172 if (output_env_get_print_xvgr_codes(oenv))
1174 fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", gn_fit);
1176 xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
1178 for (i = 0; (i < teller); i++)
1180 if (bSplit && i > 0 && abs(time[i]) < 1e-5)
1184 fprintf(fp, "%12.7f", time[i]);
1185 for (j = 0; (j < nrms); j++)
1187 fprintf(fp, " %12.7f", rlsm[j][i]);
1196 sprintf(buf, "Average %s", whatxvgname[ewhat]);
1197 sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
1198 fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
1199 for (j = 0; (j < nrms); j++)
1201 fprintf(fp, "%10d %10g\n", j, rlstot/teller);
1208 fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
1209 for (j = 0; (j < irms[0]); j++)
1211 fprintf(fp, "%10d %10g\n", j, rlsnorm[j]/teller);
1215 do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
1216 do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);
1217 do_view(oenv, opt2fn_null("-mir", NFILE, fnm), NULL);
1218 do_view(oenv, opt2fn_null("-m", NFILE, fnm), NULL);
1219 do_view(oenv, opt2fn_null("-bm", NFILE, fnm), NULL);
1220 do_view(oenv, opt2fn_null("-dist", NFILE, fnm), NULL);