3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
49 #include "gmx_fatal.h"
61 static void norm_princ(t_atoms *atoms, int isize, atom_id *index, int natoms,
67 /* equalize principal components: */
68 /* orient principal axes, get principal components */
69 orient_princ(atoms, isize, index, natoms, x, NULL, princ);
71 /* calc our own principal components */
73 for (m = 0; m < DIM; m++)
75 for (i = 0; i < isize; i++)
76 vec[m] += sqr(x[index[i]][m]);
77 vec[m] = sqrt(vec[m] / isize);
78 /* calculate scaling constants */
79 vec[m] = 1 / (sqrt(3) * vec[m]);
82 /* scale coordinates */
83 for (i = 0; i < natoms; i++)
85 for (m = 0; m < DIM; m++)
92 int gmx_rms(int argc, char *argv[])
97 "[TT]g_rms[tt] compares two structures by computing the root mean square",
98 "deviation (RMSD), the size-independent 'rho' similarity parameter",
99 "(rho) or the scaled rho (rhosc), ",
100 "see Maiorov & Crippen, PROTEINS [BB]22[bb], 273 (1995).",
101 "This is selected by [TT]-what[tt].[PAR]"
103 "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
104 "reference structure. The reference structure",
105 "is taken from the structure file ([TT]-s[tt]).[PAR]",
107 "With option [TT]-mir[tt] also a comparison with the mirror image of",
108 "the reference structure is calculated.",
109 "This is useful as a reference for 'significant' values, see",
110 "Maiorov & Crippen, PROTEINS [BB]22[bb], 273 (1995).[PAR]",
112 "Option [TT]-prev[tt] produces the comparison with a previous frame",
113 "the specified number of frames ago.[PAR]",
115 "Option [TT]-m[tt] produces a matrix in [TT].xpm[tt] format of",
116 "comparison values of each structure in the trajectory with respect to",
117 "each other structure. This file can be visualized with for instance",
118 "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].[PAR]",
120 "Option [TT]-fit[tt] controls the least-squares fitting of",
121 "the structures on top of each other: complete fit (rotation and",
122 "translation), translation only, or no fitting at all.[PAR]",
124 "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
125 "If you select the option (default) and ",
126 "supply a valid [TT].tpr[tt] file masses will be taken from there, ",
127 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
128 "[TT]GMXLIB[tt]. This is fine for proteins, but not",
129 "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
130 "is assigned to unknown atoms. You can check whether this happend by",
131 "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
133 "With [TT]-f2[tt], the 'other structures' are taken from a second",
134 "trajectory, this generates a comparison matrix of one trajectory",
135 "versus the other.[PAR]",
137 "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
139 "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
140 "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
141 "comparison group are considered." };
142 static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
143 static gmx_bool bDeltaLog = FALSE;
144 static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
145 static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1,
146 bond_user_min = -1, delta_maxy = 0.0;
147 /* strings and things for selecting difference method */
150 ewSel, ewRMSD, ewRho, ewRhoSc, ewNR
153 const char *what[ewNR + 1] =
154 { NULL, "rmsd", "rho", "rhosc", NULL };
155 const char *whatname[ewNR] =
156 { NULL, "RMSD", "Rho", "Rho sc" };
157 const char *whatlabel[ewNR] =
158 { NULL, "RMSD (nm)", "Rho", "Rho sc" };
159 const char *whatxvgname[ewNR] =
160 { NULL, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
161 const char *whatxvglabel[ewNR] =
162 { NULL, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
163 /* strings and things for fitting methods */
166 efSel, efFit, efReset, efNone, efNR
169 const char *fit[efNR + 1] =
170 { NULL, "rot+trans", "translation", "none", NULL };
171 const char *fitgraphlabel[efNR + 1] =
172 { NULL, "lsq fit", "translational fit", "no fit" };
174 static gmx_bool bMassWeighted = TRUE;
177 { "-what", FALSE, etENUM,
178 { what }, "Structural difference measure" },
179 { "-pbc", FALSE, etBOOL,
180 { &bPBC }, "PBC check" },
181 { "-fit", FALSE, etENUM,
182 { fit }, "Fit to reference structure" },
183 { "-prev", FALSE, etINT,
184 { &prev }, "Compare with previous frame" },
185 { "-split", FALSE, etBOOL,
186 { &bSplit }, "Split graph where time is zero" },
187 { "-fitall", FALSE, etBOOL,
188 { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
189 { "-skip", FALSE, etINT,
190 { &freq }, "Only write every nr-th frame to matrix" },
191 { "-skip2", FALSE, etINT,
192 { &freq2 }, "Only write every nr-th frame to matrix" },
193 { "-max", FALSE, etREAL,
194 { &rmsd_user_max }, "Maximum level in comparison matrix" },
195 { "-min", FALSE, etREAL,
196 { &rmsd_user_min }, "Minimum level in comparison matrix" },
197 { "-bmax", FALSE, etREAL,
198 { &bond_user_max }, "Maximum level in bond angle matrix" },
199 { "-bmin", FALSE, etREAL,
200 { &bond_user_min }, "Minimum level in bond angle matrix" },
201 { "-mw", FALSE, etBOOL,
202 { &bMassWeighted }, "Use mass weighting for superposition" },
203 { "-nlevels", FALSE, etINT,
204 { &nlevels }, "Number of levels in the matrices" },
205 { "-ng", FALSE, etINT,
206 { &nrms }, "Number of groups to compute RMS between" },
207 { "-dlog", FALSE, etBOOL,
209 "HIDDENUse a log x-axis in the delta t matrix" },
210 { "-dmax", FALSE, etREAL,
211 { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
212 { "-aver", FALSE, etINT,
214 "HIDDENAverage over this distance in the RMSD matrix" } };
215 int natoms_trx, natoms_trx2, natoms;
216 int i, j, k, m, teller, teller2, tel_mat, tel_mat2;
218 int maxframe = NFRAME, maxframe2 = NFRAME;
219 real t, *w_rls, *w_rms, *w_rls_m = NULL, *w_rms_m = NULL;
220 gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
221 gmx_bool bFit, bReset;
224 t_iatom *iatom = NULL;
227 rvec *x, *xp, *xm = NULL, **mat_x = NULL, **mat_x2, *mat_x2_j = NULL, vec1,
230 char buf[256], buf2[256];
233 real rlstot = 0, **rls, **rlsm = NULL, *time, *time2, *rlsnorm = NULL,
234 **rmsd_mat = NULL, **bond_mat = NULL, *axis, *axis2, *del_xaxis,
235 *del_yaxis, rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
236 real **rmsdav_mat = NULL, av_tot, weight, weight_tot;
237 real **delta = NULL, delta_max, delta_scalex = 0, delta_scaley = 0,
239 int delta_xsize = 0, del_lev = 100, mx, my, abs_my;
240 gmx_bool bA1, bA2, bPrev, bTop, *bInMat = NULL;
241 int ifit, *irms, ibond = 0, *ind_bond1 = NULL, *ind_bond2 = NULL, n_ind_m =
243 atom_id *ind_fit, **ind_rms, *ind_m = NULL, *rev_ind_m = NULL, *ind_rms_m =
245 char *gn_fit, **gn_rms;
248 gmx_rmpbc_t gpbc=NULL;
252 { efTPS, NULL, NULL, ffREAD },
253 { efTRX, "-f", NULL, ffREAD },
254 { efTRX, "-f2", NULL, ffOPTRD },
255 { efNDX, NULL, NULL, ffOPTRD },
256 { efXVG, NULL, "rmsd", ffWRITE },
257 { efXVG, "-mir", "rmsdmir", ffOPTWR },
258 { efXVG, "-a", "avgrp", ffOPTWR },
259 { efXVG, "-dist", "rmsd-dist", ffOPTWR },
260 { efXPM, "-m", "rmsd", ffOPTWR },
261 { efDAT, "-bin", "rmsd", ffOPTWR },
262 { efXPM, "-bm", "bond", ffOPTWR } };
263 #define NFILE asize(fnm)
265 CopyRight(stderr, argv[0]);
266 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW
267 | PCA_BE_NICE, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
269 /* parse enumerated options: */
271 if (ewhat == ewRho || ewhat == ewRhoSc)
272 please_cite(stdout, "Maiorov95");
274 bFit = efit == efFit;
275 bReset = efit == efReset;
278 bReset = TRUE; /* for fit, reset *must* be set */
285 /* mark active cmdline options */
286 bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
287 bFile2 = opt2bSet("-f2", NFILE, fnm);
288 bMat = opt2bSet("-m", NFILE, fnm);
289 bBond = opt2bSet("-bm", NFILE, fnm);
290 bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
291 * your RMSD matrix (hidden option */
292 bNorm = opt2bSet("-a", NFILE, fnm);
293 bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
296 fprintf(stderr, "The number of frames to skip is <= 0. "
297 "Writing out all frames.\n\n");
304 else if (bFile2 && freq2 <= 0)
307 "The number of frames to skip in second trajectory is <= 0.\n"
308 " Writing out all frames.\n\n");
317 fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
320 if (bFile2 && !bMat && !bBond)
324 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
325 " will not read from %s\n", opt2fn("-f2", NFILE,
336 "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
337 " will not read from %s\n", opt2fn("-f2",
343 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xp,
345 snew(w_rls,top.atoms.nr);
346 snew(w_rms,top.atoms.nr);
351 "WARNING: Need a run input file for bond angle matrix,\n"
352 " will not calculate bond angle matrix.\n");
358 fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares"
360 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit,
368 if (bFit && ifit < 3)
369 gmx_fatal(FARGS,"Need >= 3 points to fit!\n" );
372 for(i=0; i<ifit; i++)
375 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
377 w_rls[ind_fit[i]] = 1;
378 bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
382 fprintf(stderr,"All masses in the fit group are 0, using masses of 1\n");
383 for(i=0; i<ifit; i++)
385 w_rls[ind_fit[i]] = 1;
397 fprintf(stderr,"Select group%s for %s calculation\n",
398 (nrms>1) ? "s" : "",whatname[ewhat]);
399 get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
400 nrms,irms,ind_rms,gn_rms);
404 snew(rlsnorm,irms[0]);
407 for(j=0; j<nrms; j++)
408 snew(rls[j],maxframe);
412 for(j=0; j<nrms; j++)
413 snew(rlsm[j],maxframe);
416 for(j=0; j<nrms; j++)
419 for(i=0; i<irms[j]; i++)
422 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
424 w_rms[ind_rms[j][i]] = 1.0;
425 bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
428 fprintf(stderr,"All masses in group %d are 0, using masses of 1\n",j);
429 for(i=0; i<irms[j]; i++)
430 w_rms[ind_rms[j][i]] = 1;
433 /* Prepare reference frame */
435 gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
436 gmx_rmpbc(gpbc,top.atoms.nr,box,xp);
439 reset_x(ifit,ind_fit,top.atoms.nr,NULL,xp,w_rls);
441 /* generate reference structure mirror image: */
442 snew(xm, top.atoms.nr);
443 for(i=0; i<top.atoms.nr; i++) {
444 copy_rvec(xp[i],xm[i]);
445 xm[i][XX] = -xm[i][XX];
449 norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
451 /* read first frame */
452 natoms_trx=read_first_x(oenv,&status,opt2fn("-f",NFILE,fnm),&t,&x,box);
453 if (natoms_trx != top.atoms.nr)
455 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
456 top.atoms.nr,natoms_trx);
457 natoms = min(top.atoms.nr,natoms_trx);
458 if (bMat || bBond || bPrev) {
462 /* With -prev we use all atoms for simplicity */
465 /* Check which atoms we need (fit/rms) */
467 for(i=0; i<ifit; i++)
468 bInMat[ind_fit[i]] = TRUE;
470 for(i=0; i<irms[0]; i++)
471 if (!bInMat[ind_rms[0][i]]) {
472 bInMat[ind_rms[0][i]] = TRUE;
476 /* Make an index of needed atoms */
478 snew(rev_ind_m,natoms);
480 for(i=0; i<natoms; i++)
481 if (bPrev || bInMat[i]) {
486 snew(w_rls_m,n_ind_m);
487 snew(ind_rms_m,irms[0]);
488 snew(w_rms_m,n_ind_m);
489 for(i=0; i<ifit; i++)
490 w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
491 for(i=0; i<irms[0]; i++) {
492 ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
493 w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
500 for(k=0; k<F_NRE; k++)
501 if (IS_CHEMBOND(k)) {
502 iatom = top.idef.il[k].iatoms;
503 ncons += top.idef.il[k].nr/3;
505 fprintf(stderr,"Found %d bonds in topology\n",ncons);
506 snew(ind_bond1,ncons);
507 snew(ind_bond2,ncons);
509 for(k=0; k<F_NRE; k++)
510 if (IS_CHEMBOND(k)) {
511 iatom = top.idef.il[k].iatoms;
512 ncons = top.idef.il[k].nr/3;
513 for (i=0; i<ncons; i++) {
516 for (j=0; j<irms[0]; j++) {
517 if (iatom[3*i+1] == ind_rms[0][j]) bA1=TRUE;
518 if (iatom[3*i+2] == ind_rms[0][j]) bA2=TRUE;
521 ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
522 ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
527 fprintf(stderr,"Using %d bonds for bond angle matrix\n",ibond);
529 gmx_fatal(FARGS,"0 bonds found");
532 /* start looping over frames: */
537 gmx_rmpbc(gpbc,natoms,box,x);
540 reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
542 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
545 /*do the least squares fit to original structure*/
546 do_fit(natoms,w_rls,xp,x);
548 if (teller % freq == 0) {
549 /* keep frame for matrix calculation */
550 if (bMat || bBond || bPrev) {
551 if (tel_mat >= NFRAME)
552 srenew(mat_x,tel_mat+1);
553 snew(mat_x[tel_mat],n_ind_m);
554 for (i=0;i<n_ind_m;i++)
555 copy_rvec(x[ind_m[i]],mat_x[tel_mat][i]);
560 /*calculate energy of root_least_squares*/
565 for (i=0;i<n_ind_m;i++)
566 copy_rvec(mat_x[j][i],xp[ind_m[i]]);
568 reset_x(ifit,ind_fit,natoms,NULL,xp,w_rls);
570 do_fit(natoms,w_rls,x,xp);
572 for(j=0; (j<nrms); j++)
574 calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xp);
576 for(j=0; (j<irms[0]); j++)
578 calc_similar_ind(ewhat!=ewRMSD,1,&(ind_rms[0][j]),w_rms,x,xp);
583 /*do the least squares fit to mirror of original structure*/
584 do_fit(natoms,w_rls,xm,x);
586 for(j=0; j<nrms; j++)
588 calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xm);
590 time[teller]=output_env_conv_time(oenv,t);
593 if (teller >= maxframe) {
595 srenew(time,maxframe);
596 for(j=0; (j<nrms); j++)
597 srenew(rls[j],maxframe);
599 for(j=0; (j<nrms); j++)
600 srenew(rlsm[j],maxframe);
602 } while (read_next_x(oenv,status,&t,natoms_trx,x,box));
606 snew(time2,maxframe2);
608 fprintf(stderr,"\nWill read second trajectory file\n");
611 read_first_x(oenv,&status,opt2fn("-f2",NFILE,fnm),&t,&x,box);
612 if ( natoms_trx2 != natoms_trx )
614 "Second trajectory (%d atoms) does not match the first one"
615 " (%d atoms)", natoms_trx2, natoms_trx);
620 gmx_rmpbc(gpbc,natoms,box,x);
623 reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
625 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);
631 if (teller2 % freq2 == 0) {
632 /* keep frame for matrix calculation */
634 if (tel_mat2 >= NFRAME)
635 srenew(mat_x2,tel_mat2+1);
636 snew(mat_x2[tel_mat2],n_ind_m);
637 for (i=0;i<n_ind_m;i++)
638 copy_rvec(x[ind_m[i]],mat_x2[tel_mat2][i]);
643 time2[teller2]=output_env_conv_time(oenv,t);
646 if (teller2 >= maxframe2) {
648 srenew(time2,maxframe2);
650 } while (read_next_x(oenv,status,&t,natoms_trx2,x,box));
658 gmx_rmpbc_done(gpbc);
661 /* calculate RMS matrix */
662 fprintf(stderr,"\n");
664 fprintf(stderr,"Building %s matrix, %dx%d elements\n",
665 whatname[ewhat],tel_mat,tel_mat2);
666 snew(rmsd_mat,tel_mat);
669 fprintf(stderr,"Building bond angle matrix, %dx%d elements\n",
671 snew(bond_mat,tel_mat);
674 snew(axis2,tel_mat2);
683 for(j=0; j<tel_mat2; j++)
684 axis2[j]=time2[freq2*j];
687 delta_scalex=8.0/log(2.0);
688 delta_xsize=(int)(log(tel_mat/2)*delta_scalex+0.5)+1;
691 delta_xsize=tel_mat/2;
693 delta_scaley=1.0/delta_maxy;
694 snew(delta,delta_xsize);
695 for(j=0; j<delta_xsize; j++)
696 snew(delta[j],del_lev+1);
698 snew(rmsdav_mat,tel_mat);
699 for(j=0; j<tel_mat; j++)
700 snew(rmsdav_mat[j],tel_mat);
705 snew(mat_x2_j,natoms);
706 for(i=0; i<tel_mat; i++) {
707 axis[i]=time[freq*i];
708 fprintf(stderr,"\r element %5d; time %5.2f ",i,axis[i]);
709 if (bMat) snew(rmsd_mat[i],tel_mat2);
710 if (bBond) snew(bond_mat[i],tel_mat2);
711 for(j=0; j<tel_mat2; j++) {
713 for (k=0;k<n_ind_m;k++)
714 copy_rvec(mat_x2[j][k],mat_x2_j[k]);
715 do_fit(n_ind_m,w_rls_m,mat_x[i],mat_x2_j);
719 if (bFile2 || (i<j)) {
721 calc_similar_ind(ewhat!=ewRMSD,irms[0],ind_rms_m,
722 w_rms_m,mat_x[i],mat_x2_j);
723 if (rmsd_mat[i][j] > rmsd_max) rmsd_max=rmsd_mat[i][j];
724 if (rmsd_mat[i][j] < rmsd_min) rmsd_min=rmsd_mat[i][j];
725 rmsd_avg += rmsd_mat[i][j];
727 rmsd_mat[i][j]=rmsd_mat[j][i];
730 if (bFile2 || (i<=j)) {
732 for(m=0;m<ibond;m++) {
733 rvec_sub(mat_x[i][ind_bond1[m]],mat_x[i][ind_bond2[m]],vec1);
734 rvec_sub(mat_x2_j[ind_bond1[m]],mat_x2_j[ind_bond2[m]],vec2);
735 ang += acos(cos_angle(vec1,vec2));
737 bond_mat[i][j]=ang*180.0/(M_PI*ibond);
738 if (bond_mat[i][j] > bond_max) bond_max=bond_mat[i][j];
739 if (bond_mat[i][j] < bond_min) bond_min=bond_mat[i][j];
742 bond_mat[i][j]=bond_mat[j][i];
747 rmsd_avg /= tel_mat*tel_mat2;
749 rmsd_avg /= tel_mat*(tel_mat - 1)/2;
750 if (bMat && (avl > 0)) {
754 for(j=0; j<tel_mat-1; j++) {
755 for(i=j+1; i<tel_mat; i++) {
758 for (my=-avl; my<=avl; my++) {
759 if ((j+my>=0) && (j+my<tel_mat)) {
761 for (mx=-avl; mx<=avl; mx++) {
762 if ((i+mx>=0) && (i+mx<tel_mat)) {
763 weight = (real)(avl+1-max(abs(mx),abs_my));
764 av_tot += weight*rmsd_mat[i+mx][j+my];
770 rmsdav_mat[i][j] = av_tot/weight_tot;
771 rmsdav_mat[j][i] = rmsdav_mat[i][j];
772 if (rmsdav_mat[i][j] > rmsd_max) rmsd_max=rmsdav_mat[i][j];
779 fprintf(stderr,"\n%s: Min %f, Max %f, Avg %f\n",
780 whatname[ewhat],rmsd_min,rmsd_max,rmsd_avg);
781 rlo.r = 1; rlo.g = 1; rlo.b = 1;
782 rhi.r = 0; rhi.g = 0; rhi.b = 0;
783 if (rmsd_user_max != -1) rmsd_max=rmsd_user_max;
784 if (rmsd_user_min != -1) rmsd_min=rmsd_user_min;
785 if ((rmsd_user_max != -1) || (rmsd_user_min != -1))
786 fprintf(stderr,"Min and Max value set to resp. %f and %f\n",
788 sprintf(buf,"%s %s matrix",gn_rms[0],whatname[ewhat]);
789 write_xpm(opt2FILE("-m",NFILE,fnm,"w"),0,buf,whatlabel[ewhat],
790 output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
791 axis,axis2, rmsd_mat,rmsd_min,rmsd_max,rlo,rhi,&nlevels);
792 /* Print the distribution of RMSD values */
793 if (opt2bSet("-dist",NFILE,fnm))
794 low_rmsd_dist(opt2fn("-dist",NFILE,fnm),rmsd_max,tel_mat,rmsd_mat,oenv);
797 snew(delta_tot,delta_xsize);
798 for(j=0; j<tel_mat-1; j++) {
799 for(i=j+1; i<tel_mat; i++) {
801 if (mx < tel_mat/2) {
803 mx=(int)(log(mx)*delta_scalex+0.5);
804 my=(int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
805 delta_tot[mx] += 1.0;
806 if ((rmsd_mat[i][j]>=0) && (rmsd_mat[i][j]<=delta_maxy))
807 delta[mx][my] += 1.0;
812 for(i=0; i<delta_xsize; i++) {
813 if (delta_tot[i] > 0.0) {
814 delta_tot[i]=1.0/delta_tot[i];
815 for(j=0; j<=del_lev; j++) {
816 delta[i][j] *= delta_tot[i];
817 if (delta[i][j] > delta_max)
818 delta_max=delta[i][j];
822 fprintf(stderr,"Maximum in delta matrix: %f\n",delta_max);
823 snew(del_xaxis,delta_xsize);
824 snew(del_yaxis,del_lev+1);
825 for (i=0; i<delta_xsize; i++)
826 del_xaxis[i]=axis[i]-axis[0];
827 for (i=0; i<del_lev+1; i++)
828 del_yaxis[i]=delta_maxy*i/del_lev;
829 sprintf(buf,"%s %s vs. delta t",gn_rms[0],whatname[ewhat]);
830 fp = ffopen("delta.xpm","w");
831 write_xpm(fp,0,buf,"density",output_env_get_time_label(oenv),whatlabel[ewhat],
832 delta_xsize,del_lev+1,del_xaxis,del_yaxis,
833 delta,0.0,delta_max,rlo,rhi,&nlevels);
836 if (opt2bSet("-bin",NFILE,fnm)) {
837 /* NB: File must be binary if we use fwrite */
838 fp=ftp2FILE(efDAT,NFILE,fnm,"wb");
839 for(i=0;i<tel_mat;i++)
840 if(fwrite(rmsd_mat[i],sizeof(**rmsd_mat),tel_mat2,fp) != tel_mat2)
842 gmx_fatal(FARGS,"Error writing to output file");
848 fprintf(stderr,"\nMin. angle: %f, Max. angle: %f\n",bond_min,bond_max);
849 if (bond_user_max != -1) bond_max=bond_user_max;
850 if (bond_user_min != -1) bond_min=bond_user_min;
851 if ((bond_user_max != -1) || (bond_user_min != -1))
852 fprintf(stderr,"Bond angle Min and Max set to:\n"
853 "Min. angle: %f, Max. angle: %f\n",bond_min,bond_max);
854 rlo.r = 1; rlo.g = 1; rlo.b = 1;
855 rhi.r = 0; rhi.g = 0; rhi.b = 0;
856 sprintf(buf,"%s av. bond angle deviation",gn_rms[0]);
857 write_xpm(opt2FILE("-bm",NFILE,fnm,"w"),0,buf,"degrees",
858 output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
859 axis,axis2, bond_mat,bond_min,bond_max,rlo,rhi,&nlevels);
863 bAv=opt2bSet("-a",NFILE,fnm);
865 /* Write the RMSD's to file */
867 sprintf(buf,"%s",whatxvgname[ewhat]);
869 sprintf(buf,"%s with frame %g %s ago",whatxvgname[ewhat],
870 time[prev*freq]-time[0], output_env_get_time_label(oenv));
871 fp=xvgropen(opt2fn("-o",NFILE,fnm),buf,output_env_get_xvgr_tlabel(oenv),
872 whatxvglabel[ewhat], oenv);
873 if (output_env_get_print_xvgr_codes(oenv))
874 fprintf(fp,"@ subtitle \"%s%s after %s%s%s\"\n",
875 (nrms==1)?"":"of " , gn_rms[0], fitgraphlabel[efit],
876 bFit ?" to ":"" , bFit?gn_fit:"");
878 xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
879 for(i=0; (i<teller); i++) {
880 if ( bSplit && i>0 &&
881 abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv))<1e-5 )
885 fprintf(fp,"%12.7f",time[bPrev ? freq*i : i]);
886 for(j=0; (j<nrms); j++) {
887 fprintf(fp," %12.7f",rls[j][i]);
896 /* Write the mirror RMSD's to file */
897 sprintf(buf,"%s with Mirror",whatxvgname[ewhat]);
898 sprintf(buf2,"Mirror %s",whatxvglabel[ewhat]);
899 fp=xvgropen(opt2fn("-mir",NFILE,fnm), buf, output_env_get_xvgr_tlabel(oenv),
902 if (output_env_get_print_xvgr_codes(oenv))
903 fprintf(fp,"@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
907 if (output_env_get_print_xvgr_codes(oenv))
908 fprintf(fp,"@ subtitle \"after lsq fit to mirror %s\"\n",gn_fit);
909 xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
911 for(i=0; (i<teller); i++) {
912 if ( bSplit && i>0 && abs(time[i])<1e-5 )
914 fprintf(fp,"%12.7f",time[i]);
915 for(j=0; (j<nrms); j++)
916 fprintf(fp," %12.7f",rlsm[j][i]);
923 sprintf(buf,"Average %s",whatxvgname[ewhat]);
924 sprintf(buf2,"Average %s",whatxvglabel[ewhat]);
925 fp = xvgropen(opt2fn("-a",NFILE,fnm), buf, "Residue", buf2,oenv);
926 for(j=0; (j<nrms); j++)
927 fprintf(fp,"%10d %10g\n",j,rlstot/teller);
932 fp = xvgropen("aver.xvg",gn_rms[0],"Residue",whatxvglabel[ewhat],oenv);
933 for(j=0; (j<irms[0]); j++)
934 fprintf(fp,"%10d %10g\n",j,rlsnorm[j]/teller);
937 do_view(oenv,opt2fn_null("-a",NFILE,fnm),"-graphtype bar");
938 do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
939 do_view(oenv,opt2fn_null("-mir",NFILE,fnm),NULL);
940 do_view(oenv,opt2fn_null("-m",NFILE,fnm),NULL);
941 do_view(oenv,opt2fn_null("-bm",NFILE,fnm),NULL);
942 do_view(oenv,opt2fn_null("-dist",NFILE,fnm),NULL);