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,2015,2016, 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.
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/matio.h"
44 #include "gromacs/fileio/tpxio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/cmat.h"
48 #include "gromacs/gmxana/gmx_ana.h"
49 #include "gromacs/gmxana/princ.h"
50 #include "gromacs/legacyheaders/copyrite.h"
51 #include "gromacs/legacyheaders/macros.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/legacyheaders/viewit.h"
54 #include "gromacs/math/do_fit.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/pbcutil/rmpbc.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
62 static void norm_princ(t_atoms
*atoms
, int isize
, atom_id
*index
, int natoms
,
68 /* equalize principal components: */
69 /* orient principal axes, get principal components */
70 orient_princ(atoms
, isize
, index
, natoms
, x
, NULL
, princ
);
72 /* calc our own principal components */
74 for (m
= 0; m
< DIM
; m
++)
76 for (i
= 0; i
< isize
; i
++)
78 vec
[m
] += sqr(x
[index
[i
]][m
]);
80 vec
[m
] = sqrt(vec
[m
] / isize
);
81 /* calculate scaling constants */
82 vec
[m
] = 1 / (sqrt(3) * vec
[m
]);
85 /* scale coordinates */
86 for (i
= 0; i
< natoms
; i
++)
88 for (m
= 0; m
< DIM
; m
++)
95 int gmx_rms(int argc
, char *argv
[])
99 "[THISMODULE] compares two structures by computing the root mean square",
100 "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
101 "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
102 "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
103 "This is selected by [TT]-what[tt].[PAR]"
105 "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
106 "reference structure. The reference structure",
107 "is taken from the structure file ([TT]-s[tt]).[PAR]",
109 "With option [TT]-mir[tt] also a comparison with the mirror image of",
110 "the reference structure is calculated.",
111 "This is useful as a reference for 'significant' values, see",
112 "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
114 "Option [TT]-prev[tt] produces the comparison with a previous frame",
115 "the specified number of frames ago.[PAR]",
117 "Option [TT]-m[tt] produces a matrix in [REF].xpm[ref] format of",
118 "comparison values of each structure in the trajectory with respect to",
119 "each other structure. This file can be visualized with for instance",
120 "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]",
122 "Option [TT]-fit[tt] controls the least-squares fitting of",
123 "the structures on top of each other: complete fit (rotation and",
124 "translation), translation only, or no fitting at all.[PAR]",
126 "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
127 "If you select the option (default) and ",
128 "supply a valid [REF].tpr[ref] file masses will be taken from there, ",
129 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
130 "[TT]GMXLIB[tt]. This is fine for proteins, but not",
131 "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
132 "is assigned to unknown atoms. You can check whether this happened by",
133 "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
135 "With [TT]-f2[tt], the 'other structures' are taken from a second",
136 "trajectory, this generates a comparison matrix of one trajectory",
137 "versus the other.[PAR]",
139 "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
141 "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
142 "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
143 "comparison group are considered."
145 static gmx_bool bPBC
= TRUE
, bFitAll
= TRUE
, bSplit
= FALSE
;
146 static gmx_bool bDeltaLog
= FALSE
;
147 static int prev
= 0, freq
= 1, freq2
= 1, nlevels
= 80, avl
= 0;
148 static real rmsd_user_max
= -1, rmsd_user_min
= -1, bond_user_max
= -1,
149 bond_user_min
= -1, delta_maxy
= 0.0;
150 /* strings and things for selecting difference method */
153 ewSel
, ewRMSD
, ewRho
, ewRhoSc
, ewNR
156 const char *what
[ewNR
+ 1] =
157 { NULL
, "rmsd", "rho", "rhosc", NULL
};
158 const char *whatname
[ewNR
] =
159 { NULL
, "RMSD", "Rho", "Rho sc" };
160 const char *whatlabel
[ewNR
] =
161 { NULL
, "RMSD (nm)", "Rho", "Rho sc" };
162 const char *whatxvgname
[ewNR
] =
163 { NULL
, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
164 const char *whatxvglabel
[ewNR
] =
165 { NULL
, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
166 /* strings and things for fitting methods */
169 efSel
, efFit
, efReset
, efNone
, efNR
172 const char *fit
[efNR
+ 1] =
173 { NULL
, "rot+trans", "translation", "none", NULL
};
174 const char *fitgraphlabel
[efNR
+ 1] =
175 { NULL
, "lsq fit", "translational fit", "no fit" };
177 static gmx_bool bMassWeighted
= TRUE
;
180 { "-what", FALSE
, etENUM
,
181 { what
}, "Structural difference measure" },
182 { "-pbc", FALSE
, etBOOL
,
183 { &bPBC
}, "PBC check" },
184 { "-fit", FALSE
, etENUM
,
185 { fit
}, "Fit to reference structure" },
186 { "-prev", FALSE
, etINT
,
187 { &prev
}, "Compare with previous frame" },
188 { "-split", FALSE
, etBOOL
,
189 { &bSplit
}, "Split graph where time is zero" },
190 { "-fitall", FALSE
, etBOOL
,
191 { &bFitAll
}, "HIDDENFit all pairs of structures in matrix" },
192 { "-skip", FALSE
, etINT
,
193 { &freq
}, "Only write every nr-th frame to matrix" },
194 { "-skip2", FALSE
, etINT
,
195 { &freq2
}, "Only write every nr-th frame to matrix" },
196 { "-max", FALSE
, etREAL
,
197 { &rmsd_user_max
}, "Maximum level in comparison matrix" },
198 { "-min", FALSE
, etREAL
,
199 { &rmsd_user_min
}, "Minimum level in comparison matrix" },
200 { "-bmax", FALSE
, etREAL
,
201 { &bond_user_max
}, "Maximum level in bond angle matrix" },
202 { "-bmin", FALSE
, etREAL
,
203 { &bond_user_min
}, "Minimum level in bond angle matrix" },
204 { "-mw", FALSE
, etBOOL
,
205 { &bMassWeighted
}, "Use mass weighting for superposition" },
206 { "-nlevels", FALSE
, etINT
,
207 { &nlevels
}, "Number of levels in the matrices" },
208 { "-ng", FALSE
, etINT
,
209 { &nrms
}, "Number of groups to compute RMS between" },
210 { "-dlog", FALSE
, etBOOL
,
212 "HIDDENUse a log x-axis in the delta t matrix" },
213 { "-dmax", FALSE
, etREAL
,
214 { &delta_maxy
}, "HIDDENMaximum level in delta matrix" },
215 { "-aver", FALSE
, etINT
,
217 "HIDDENAverage over this distance in the RMSD matrix" }
219 int natoms_trx
, natoms_trx2
, natoms
;
220 int i
, j
, k
, m
, teller
, teller2
, tel_mat
, tel_mat2
;
222 int maxframe
= NFRAME
, maxframe2
= NFRAME
;
223 real t
, *w_rls
, *w_rms
, *w_rls_m
= NULL
, *w_rms_m
= NULL
;
224 gmx_bool bNorm
, bAv
, bFreq2
, bFile2
, bMat
, bBond
, bDelta
, bMirror
, bMass
;
225 gmx_bool bFit
, bReset
;
228 t_iatom
*iatom
= NULL
;
231 rvec
*x
, *xp
, *xm
= NULL
, **mat_x
= NULL
, **mat_x2
, *mat_x2_j
= NULL
, vec1
,
234 char buf
[256], buf2
[256];
237 real rlstot
= 0, **rls
, **rlsm
= NULL
, *time
, *time2
, *rlsnorm
= NULL
,
238 **rmsd_mat
= NULL
, **bond_mat
= NULL
, *axis
, *axis2
, *del_xaxis
,
239 *del_yaxis
, rmsd_max
, rmsd_min
, rmsd_avg
, bond_max
, bond_min
, ang
;
240 real
**rmsdav_mat
= NULL
, av_tot
, weight
, weight_tot
;
241 real
**delta
= NULL
, delta_max
, delta_scalex
= 0, delta_scaley
= 0,
243 int delta_xsize
= 0, del_lev
= 100, mx
, my
, abs_my
;
244 gmx_bool bA1
, bA2
, bPrev
, bTop
, *bInMat
= NULL
;
245 int ifit
, *irms
, ibond
= 0, *ind_bond1
= NULL
, *ind_bond2
= NULL
, n_ind_m
=
247 atom_id
*ind_fit
, **ind_rms
, *ind_m
= NULL
, *rev_ind_m
= NULL
, *ind_rms_m
=
249 char *gn_fit
, **gn_rms
;
252 gmx_rmpbc_t gpbc
= NULL
;
256 { efTPS
, NULL
, NULL
, ffREAD
},
257 { efTRX
, "-f", NULL
, ffREAD
},
258 { efTRX
, "-f2", NULL
, ffOPTRD
},
259 { efNDX
, NULL
, NULL
, ffOPTRD
},
260 { efXVG
, NULL
, "rmsd", ffWRITE
},
261 { efXVG
, "-mir", "rmsdmir", ffOPTWR
},
262 { efXVG
, "-a", "avgrp", ffOPTWR
},
263 { efXVG
, "-dist", "rmsd-dist", ffOPTWR
},
264 { efXPM
, "-m", "rmsd", ffOPTWR
},
265 { efDAT
, "-bin", "rmsd", ffOPTWR
},
266 { efXPM
, "-bm", "bond", ffOPTWR
}
268 #define NFILE asize(fnm)
270 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_TIME_UNIT
| PCA_CAN_VIEW
,
271 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, 0, NULL
,
276 /* parse enumerated options: */
278 if (ewhat
== ewRho
|| ewhat
== ewRhoSc
)
280 please_cite(stdout
, "Maiorov95");
283 bFit
= efit
== efFit
;
284 bReset
= efit
== efReset
;
287 bReset
= TRUE
; /* for fit, reset *must* be set */
294 /* mark active cmdline options */
295 bMirror
= opt2bSet("-mir", NFILE
, fnm
); /* calc RMSD vs mirror of ref. */
296 bFile2
= opt2bSet("-f2", NFILE
, fnm
);
297 bMat
= opt2bSet("-m", NFILE
, fnm
);
298 bBond
= opt2bSet("-bm", NFILE
, fnm
);
299 bDelta
= (delta_maxy
> 0); /* calculate rmsd vs delta t matrix from *
300 * your RMSD matrix (hidden option */
301 bNorm
= opt2bSet("-a", NFILE
, fnm
);
302 bFreq2
= opt2parg_bSet("-skip2", asize(pa
), pa
);
305 fprintf(stderr
, "The number of frames to skip is <= 0. "
306 "Writing out all frames.\n\n");
313 else if (bFile2
&& freq2
<= 0)
316 "The number of frames to skip in second trajectory is <= 0.\n"
317 " Writing out all frames.\n\n");
324 fprintf(stderr
, "WARNING: using option -prev with large trajectories will\n"
325 " require a lot of memory and could lead to crashes\n");
329 fprintf(stderr
, "WARNING: option -skip also applies to -prev\n");
333 if (bFile2
&& !bMat
&& !bBond
)
337 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
338 " will not read from %s\n", opt2fn("-f2", NFILE
,
349 "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
350 " will not read from %s\n", opt2fn("-f2",
356 bTop
= read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), buf
, &top
, &ePBC
, &xp
,
358 snew(w_rls
, top
.atoms
.nr
);
359 snew(w_rms
, top
.atoms
.nr
);
364 "WARNING: Need a run input file for bond angle matrix,\n"
365 " will not calculate bond angle matrix.\n");
371 fprintf(stderr
, "Select group for %s fit\n", bFit
? "least squares"
373 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &ifit
,
383 if (bFit
&& ifit
< 3)
385 gmx_fatal(FARGS
, "Need >= 3 points to fit!\n" );
389 for (i
= 0; i
< ifit
; i
++)
393 w_rls
[ind_fit
[i
]] = top
.atoms
.atom
[ind_fit
[i
]].m
;
397 w_rls
[ind_fit
[i
]] = 1;
399 bMass
= bMass
|| (top
.atoms
.atom
[ind_fit
[i
]].m
!= 0);
403 fprintf(stderr
, "All masses in the fit group are 0, using masses of 1\n");
404 for (i
= 0; i
< ifit
; i
++)
406 w_rls
[ind_fit
[i
]] = 1;
420 fprintf(stderr
, "Select group%s for %s calculation\n",
421 (nrms
> 1) ? "s" : "", whatname
[ewhat
]);
422 get_index(&(top
.atoms
), ftp2fn_null(efNDX
, NFILE
, fnm
),
423 nrms
, irms
, ind_rms
, gn_rms
);
427 snew(rlsnorm
, irms
[0]);
430 for (j
= 0; j
< nrms
; j
++)
432 snew(rls
[j
], maxframe
);
437 for (j
= 0; j
< nrms
; j
++)
439 snew(rlsm
[j
], maxframe
);
442 snew(time
, maxframe
);
443 for (j
= 0; j
< nrms
; j
++)
446 for (i
= 0; i
< irms
[j
]; i
++)
450 w_rms
[ind_rms
[j
][i
]] = top
.atoms
.atom
[ind_rms
[j
][i
]].m
;
454 w_rms
[ind_rms
[j
][i
]] = 1.0;
456 bMass
= bMass
|| (top
.atoms
.atom
[ind_rms
[j
][i
]].m
!= 0);
460 fprintf(stderr
, "All masses in group %d are 0, using masses of 1\n", j
);
461 for (i
= 0; i
< irms
[j
]; i
++)
463 w_rms
[ind_rms
[j
][i
]] = 1;
467 /* Prepare reference frame */
470 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, top
.atoms
.nr
);
471 gmx_rmpbc(gpbc
, top
.atoms
.nr
, box
, xp
);
475 reset_x(ifit
, ind_fit
, top
.atoms
.nr
, NULL
, xp
, w_rls
);
479 /* generate reference structure mirror image: */
480 snew(xm
, top
.atoms
.nr
);
481 for (i
= 0; i
< top
.atoms
.nr
; i
++)
483 copy_rvec(xp
[i
], xm
[i
]);
484 xm
[i
][XX
] = -xm
[i
][XX
];
487 if (ewhat
== ewRhoSc
)
489 norm_princ(&top
.atoms
, ifit
, ind_fit
, top
.atoms
.nr
, xp
);
492 /* read first frame */
493 natoms_trx
= read_first_x(oenv
, &status
, opt2fn("-f", NFILE
, fnm
), &t
, &x
, box
);
494 if (natoms_trx
!= top
.atoms
.nr
)
497 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
498 top
.atoms
.nr
, natoms_trx
);
500 natoms
= min(top
.atoms
.nr
, natoms_trx
);
501 if (bMat
|| bBond
|| bPrev
)
507 /* With -prev we use all atoms for simplicity */
512 /* Check which atoms we need (fit/rms) */
513 snew(bInMat
, natoms
);
514 for (i
= 0; i
< ifit
; i
++)
516 bInMat
[ind_fit
[i
]] = TRUE
;
519 for (i
= 0; i
< irms
[0]; i
++)
521 if (!bInMat
[ind_rms
[0][i
]])
523 bInMat
[ind_rms
[0][i
]] = TRUE
;
528 /* Make an index of needed atoms */
529 snew(ind_m
, n_ind_m
);
530 snew(rev_ind_m
, natoms
);
532 for (i
= 0; i
< natoms
; i
++)
534 if (bPrev
|| bInMat
[i
])
541 snew(w_rls_m
, n_ind_m
);
542 snew(ind_rms_m
, irms
[0]);
543 snew(w_rms_m
, n_ind_m
);
544 for (i
= 0; i
< ifit
; i
++)
546 w_rls_m
[rev_ind_m
[ind_fit
[i
]]] = w_rls
[ind_fit
[i
]];
548 for (i
= 0; i
< irms
[0]; i
++)
550 ind_rms_m
[i
] = rev_ind_m
[ind_rms
[0][i
]];
551 w_rms_m
[ind_rms_m
[i
]] = w_rms
[ind_rms
[0][i
]];
559 for (k
= 0; k
< F_NRE
; k
++)
563 iatom
= top
.idef
.il
[k
].iatoms
;
564 ncons
+= top
.idef
.il
[k
].nr
/3;
567 fprintf(stderr
, "Found %d bonds in topology\n", ncons
);
568 snew(ind_bond1
, ncons
);
569 snew(ind_bond2
, ncons
);
571 for (k
= 0; k
< F_NRE
; k
++)
575 iatom
= top
.idef
.il
[k
].iatoms
;
576 ncons
= top
.idef
.il
[k
].nr
/3;
577 for (i
= 0; i
< ncons
; i
++)
581 for (j
= 0; j
< irms
[0]; j
++)
583 if (iatom
[3*i
+1] == ind_rms
[0][j
])
587 if (iatom
[3*i
+2] == ind_rms
[0][j
])
594 ind_bond1
[ibond
] = rev_ind_m
[iatom
[3*i
+1]];
595 ind_bond2
[ibond
] = rev_ind_m
[iatom
[3*i
+2]];
601 fprintf(stderr
, "Using %d bonds for bond angle matrix\n", ibond
);
604 gmx_fatal(FARGS
, "0 bonds found");
608 /* start looping over frames: */
615 gmx_rmpbc(gpbc
, natoms
, box
, x
);
620 reset_x(ifit
, ind_fit
, natoms
, NULL
, x
, w_rls
);
622 if (ewhat
== ewRhoSc
)
624 norm_princ(&top
.atoms
, ifit
, ind_fit
, natoms
, x
);
629 /*do the least squares fit to original structure*/
630 do_fit(natoms
, w_rls
, xp
, x
);
633 if (teller
% freq
== 0)
635 /* keep frame for matrix calculation */
636 if (bMat
|| bBond
|| bPrev
)
638 if (tel_mat
>= NFRAME
)
640 srenew(mat_x
, tel_mat
+1);
642 snew(mat_x
[tel_mat
], n_ind_m
);
643 for (i
= 0; i
< n_ind_m
; i
++)
645 copy_rvec(x
[ind_m
[i
]], mat_x
[tel_mat
][i
]);
651 /*calculate energy of root_least_squares*/
659 for (i
= 0; i
< n_ind_m
; i
++)
661 copy_rvec(mat_x
[j
][i
], xp
[ind_m
[i
]]);
665 reset_x(ifit
, ind_fit
, natoms
, NULL
, xp
, w_rls
);
669 do_fit(natoms
, w_rls
, x
, xp
);
672 for (j
= 0; (j
< nrms
); j
++)
675 calc_similar_ind(ewhat
!= ewRMSD
, irms
[j
], ind_rms
[j
], w_rms
, x
, xp
);
679 for (j
= 0; (j
< irms
[0]); j
++)
682 calc_similar_ind(ewhat
!= ewRMSD
, 1, &(ind_rms
[0][j
]), w_rms
, x
, xp
);
690 /*do the least squares fit to mirror of original structure*/
691 do_fit(natoms
, w_rls
, xm
, x
);
694 for (j
= 0; j
< nrms
; j
++)
697 calc_similar_ind(ewhat
!= ewRMSD
, irms
[j
], ind_rms
[j
], w_rms
, x
, xm
);
700 time
[teller
] = output_env_conv_time(oenv
, t
);
703 if (teller
>= maxframe
)
706 srenew(time
, maxframe
);
707 for (j
= 0; (j
< nrms
); j
++)
709 srenew(rls
[j
], maxframe
);
713 for (j
= 0; (j
< nrms
); j
++)
715 srenew(rlsm
[j
], maxframe
);
720 while (read_next_x(oenv
, status
, &t
, x
, box
));
725 snew(time2
, maxframe2
);
727 fprintf(stderr
, "\nWill read second trajectory file\n");
728 snew(mat_x2
, NFRAME
);
730 read_first_x(oenv
, &status
, opt2fn("-f2", NFILE
, fnm
), &t
, &x
, box
);
731 if (natoms_trx2
!= natoms_trx
)
734 "Second trajectory (%d atoms) does not match the first one"
735 " (%d atoms)", natoms_trx2
, natoms_trx
);
743 gmx_rmpbc(gpbc
, natoms
, box
, x
);
748 reset_x(ifit
, ind_fit
, natoms
, NULL
, x
, w_rls
);
750 if (ewhat
== ewRhoSc
)
752 norm_princ(&top
.atoms
, ifit
, ind_fit
, natoms
, x
);
757 /*do the least squares fit to original structure*/
758 do_fit(natoms
, w_rls
, xp
, x
);
761 if (teller2
% freq2
== 0)
763 /* keep frame for matrix calculation */
766 if (tel_mat2
>= NFRAME
)
768 srenew(mat_x2
, tel_mat2
+1);
770 snew(mat_x2
[tel_mat2
], n_ind_m
);
771 for (i
= 0; i
< n_ind_m
; i
++)
773 copy_rvec(x
[ind_m
[i
]], mat_x2
[tel_mat2
][i
]);
779 time2
[teller2
] = output_env_conv_time(oenv
, t
);
782 if (teller2
>= maxframe2
)
785 srenew(time2
, maxframe2
);
788 while (read_next_x(oenv
, status
, &t
, x
, box
));
798 gmx_rmpbc_done(gpbc
);
802 /* calculate RMS matrix */
803 fprintf(stderr
, "\n");
806 fprintf(stderr
, "Building %s matrix, %dx%d elements\n",
807 whatname
[ewhat
], tel_mat
, tel_mat2
);
808 snew(rmsd_mat
, tel_mat
);
812 fprintf(stderr
, "Building bond angle matrix, %dx%d elements\n",
814 snew(bond_mat
, tel_mat
);
817 snew(axis2
, tel_mat2
);
830 for (j
= 0; j
< tel_mat2
; j
++)
832 axis2
[j
] = time2
[freq2
*j
];
838 delta_scalex
= 8.0/log(2.0);
839 delta_xsize
= (int)(log(tel_mat
/2)*delta_scalex
+0.5)+1;
843 delta_xsize
= tel_mat
/2;
845 delta_scaley
= 1.0/delta_maxy
;
846 snew(delta
, delta_xsize
);
847 for (j
= 0; j
< delta_xsize
; j
++)
849 snew(delta
[j
], del_lev
+1);
853 snew(rmsdav_mat
, tel_mat
);
854 for (j
= 0; j
< tel_mat
; j
++)
856 snew(rmsdav_mat
[j
], tel_mat
);
863 snew(mat_x2_j
, natoms
);
865 for (i
= 0; i
< tel_mat
; i
++)
867 axis
[i
] = time
[freq
*i
];
868 fprintf(stderr
, "\r element %5d; time %5.2f ", i
, axis
[i
]);
871 snew(rmsd_mat
[i
], tel_mat2
);
875 snew(bond_mat
[i
], tel_mat2
);
877 for (j
= 0; j
< tel_mat2
; j
++)
881 for (k
= 0; k
< n_ind_m
; k
++)
883 copy_rvec(mat_x2
[j
][k
], mat_x2_j
[k
]);
885 do_fit(n_ind_m
, w_rls_m
, mat_x
[i
], mat_x2_j
);
889 mat_x2_j
= mat_x2
[j
];
893 if (bFile2
|| (i
< j
))
896 calc_similar_ind(ewhat
!= ewRMSD
, irms
[0], ind_rms_m
,
897 w_rms_m
, mat_x
[i
], mat_x2_j
);
898 if (rmsd_mat
[i
][j
] > rmsd_max
)
900 rmsd_max
= rmsd_mat
[i
][j
];
902 if (rmsd_mat
[i
][j
] < rmsd_min
)
904 rmsd_min
= rmsd_mat
[i
][j
];
906 rmsd_avg
+= rmsd_mat
[i
][j
];
910 rmsd_mat
[i
][j
] = rmsd_mat
[j
][i
];
915 if (bFile2
|| (i
<= j
))
918 for (m
= 0; m
< ibond
; m
++)
920 rvec_sub(mat_x
[i
][ind_bond1
[m
]], mat_x
[i
][ind_bond2
[m
]], vec1
);
921 rvec_sub(mat_x2_j
[ind_bond1
[m
]], mat_x2_j
[ind_bond2
[m
]], vec2
);
922 ang
+= acos(cos_angle(vec1
, vec2
));
924 bond_mat
[i
][j
] = ang
*180.0/(M_PI
*ibond
);
925 if (bond_mat
[i
][j
] > bond_max
)
927 bond_max
= bond_mat
[i
][j
];
929 if (bond_mat
[i
][j
] < bond_min
)
931 bond_min
= bond_mat
[i
][j
];
936 bond_mat
[i
][j
] = bond_mat
[j
][i
];
943 rmsd_avg
/= tel_mat
*tel_mat2
;
947 rmsd_avg
/= tel_mat
*(tel_mat
- 1)/2;
949 if (bMat
&& (avl
> 0))
954 for (j
= 0; j
< tel_mat
-1; j
++)
956 for (i
= j
+1; i
< tel_mat
; i
++)
960 for (my
= -avl
; my
<= avl
; my
++)
962 if ((j
+my
>= 0) && (j
+my
< tel_mat
))
965 for (mx
= -avl
; mx
<= avl
; mx
++)
967 if ((i
+mx
>= 0) && (i
+mx
< tel_mat
))
969 weight
= (real
)(avl
+1-max(abs(mx
), abs_my
));
970 av_tot
+= weight
*rmsd_mat
[i
+mx
][j
+my
];
971 weight_tot
+= weight
;
976 rmsdav_mat
[i
][j
] = av_tot
/weight_tot
;
977 rmsdav_mat
[j
][i
] = rmsdav_mat
[i
][j
];
978 if (rmsdav_mat
[i
][j
] > rmsd_max
)
980 rmsd_max
= rmsdav_mat
[i
][j
];
984 rmsd_mat
= rmsdav_mat
;
989 fprintf(stderr
, "\n%s: Min %f, Max %f, Avg %f\n",
990 whatname
[ewhat
], rmsd_min
, rmsd_max
, rmsd_avg
);
991 rlo
.r
= 1; rlo
.g
= 1; rlo
.b
= 1;
992 rhi
.r
= 0; rhi
.g
= 0; rhi
.b
= 0;
993 if (rmsd_user_max
!= -1)
995 rmsd_max
= rmsd_user_max
;
997 if (rmsd_user_min
!= -1)
999 rmsd_min
= rmsd_user_min
;
1001 if ((rmsd_user_max
!= -1) || (rmsd_user_min
!= -1))
1003 fprintf(stderr
, "Min and Max value set to resp. %f and %f\n",
1004 rmsd_min
, rmsd_max
);
1006 sprintf(buf
, "%s %s matrix", gn_rms
[0], whatname
[ewhat
]);
1007 write_xpm(opt2FILE("-m", NFILE
, fnm
, "w"), 0, buf
, whatlabel
[ewhat
],
1008 output_env_get_time_label(oenv
), output_env_get_time_label(oenv
), tel_mat
, tel_mat2
,
1009 axis
, axis2
, rmsd_mat
, rmsd_min
, rmsd_max
, rlo
, rhi
, &nlevels
);
1010 /* Print the distribution of RMSD values */
1011 if (opt2bSet("-dist", NFILE
, fnm
))
1013 low_rmsd_dist(opt2fn("-dist", NFILE
, fnm
), rmsd_max
, tel_mat
, rmsd_mat
, oenv
);
1018 snew(delta_tot
, delta_xsize
);
1019 for (j
= 0; j
< tel_mat
-1; j
++)
1021 for (i
= j
+1; i
< tel_mat
; i
++)
1028 mx
= (int)(log(mx
)*delta_scalex
+0.5);
1030 my
= (int)(rmsd_mat
[i
][j
]*delta_scaley
*del_lev
+0.5);
1031 delta_tot
[mx
] += 1.0;
1032 if ((rmsd_mat
[i
][j
] >= 0) && (rmsd_mat
[i
][j
] <= delta_maxy
))
1034 delta
[mx
][my
] += 1.0;
1040 for (i
= 0; i
< delta_xsize
; i
++)
1042 if (delta_tot
[i
] > 0.0)
1044 delta_tot
[i
] = 1.0/delta_tot
[i
];
1045 for (j
= 0; j
<= del_lev
; j
++)
1047 delta
[i
][j
] *= delta_tot
[i
];
1048 if (delta
[i
][j
] > delta_max
)
1050 delta_max
= delta
[i
][j
];
1055 fprintf(stderr
, "Maximum in delta matrix: %f\n", delta_max
);
1056 snew(del_xaxis
, delta_xsize
);
1057 snew(del_yaxis
, del_lev
+1);
1058 for (i
= 0; i
< delta_xsize
; i
++)
1060 del_xaxis
[i
] = axis
[i
]-axis
[0];
1062 for (i
= 0; i
< del_lev
+1; i
++)
1064 del_yaxis
[i
] = delta_maxy
*i
/del_lev
;
1066 sprintf(buf
, "%s %s vs. delta t", gn_rms
[0], whatname
[ewhat
]);
1067 fp
= gmx_ffopen("delta.xpm", "w");
1068 write_xpm(fp
, 0, buf
, "density", output_env_get_time_label(oenv
), whatlabel
[ewhat
],
1069 delta_xsize
, del_lev
+1, del_xaxis
, del_yaxis
,
1070 delta
, 0.0, delta_max
, rlo
, rhi
, &nlevels
);
1073 if (opt2bSet("-bin", NFILE
, fnm
))
1075 /* NB: File must be binary if we use fwrite */
1076 fp
= ftp2FILE(efDAT
, NFILE
, fnm
, "wb");
1077 for (i
= 0; i
< tel_mat
; i
++)
1079 if (fwrite(rmsd_mat
[i
], sizeof(**rmsd_mat
), tel_mat2
, fp
) != tel_mat2
)
1081 gmx_fatal(FARGS
, "Error writing to output file");
1089 fprintf(stderr
, "\nMin. angle: %f, Max. angle: %f\n", bond_min
, bond_max
);
1090 if (bond_user_max
!= -1)
1092 bond_max
= bond_user_max
;
1094 if (bond_user_min
!= -1)
1096 bond_min
= bond_user_min
;
1098 if ((bond_user_max
!= -1) || (bond_user_min
!= -1))
1100 fprintf(stderr
, "Bond angle Min and Max set to:\n"
1101 "Min. angle: %f, Max. angle: %f\n", bond_min
, bond_max
);
1103 rlo
.r
= 1; rlo
.g
= 1; rlo
.b
= 1;
1104 rhi
.r
= 0; rhi
.g
= 0; rhi
.b
= 0;
1105 sprintf(buf
, "%s av. bond angle deviation", gn_rms
[0]);
1106 write_xpm(opt2FILE("-bm", NFILE
, fnm
, "w"), 0, buf
, "degrees",
1107 output_env_get_time_label(oenv
), output_env_get_time_label(oenv
), tel_mat
, tel_mat2
,
1108 axis
, axis2
, bond_mat
, bond_min
, bond_max
, rlo
, rhi
, &nlevels
);
1112 bAv
= opt2bSet("-a", NFILE
, fnm
);
1114 /* Write the RMSD's to file */
1117 sprintf(buf
, "%s", whatxvgname
[ewhat
]);
1121 sprintf(buf
, "%s with frame %g %s ago", whatxvgname
[ewhat
],
1122 time
[prev
*freq
]-time
[0], output_env_get_time_label(oenv
));
1124 fp
= xvgropen(opt2fn("-o", NFILE
, fnm
), buf
, output_env_get_xvgr_tlabel(oenv
),
1125 whatxvglabel
[ewhat
], oenv
);
1126 if (output_env_get_print_xvgr_codes(oenv
))
1128 fprintf(fp
, "@ subtitle \"%s%s after %s%s%s\"\n",
1129 (nrms
== 1) ? "" : "of ", gn_rms
[0], fitgraphlabel
[efit
],
1130 bFit
? " to " : "", bFit
? gn_fit
: "");
1134 xvgr_legend(fp
, nrms
, (const char**)gn_rms
, oenv
);
1136 for (i
= 0; (i
< teller
); i
++)
1138 if (bSplit
&& i
> 0 &&
1139 fabs(time
[bPrev
? freq
*i
: i
]/output_env_get_time_factor(oenv
)) < 1e-5)
1141 fprintf(fp
, "%s\n", output_env_get_print_xvgr_codes(oenv
) ? "&" : "");
1143 fprintf(fp
, "%12.7f", time
[bPrev
? freq
*i
: i
]);
1144 for (j
= 0; (j
< nrms
); j
++)
1146 fprintf(fp
, " %12.7f", rls
[j
][i
]);
1149 rlstot
+= rls
[j
][i
];
1158 /* Write the mirror RMSD's to file */
1159 sprintf(buf
, "%s with Mirror", whatxvgname
[ewhat
]);
1160 sprintf(buf2
, "Mirror %s", whatxvglabel
[ewhat
]);
1161 fp
= xvgropen(opt2fn("-mir", NFILE
, fnm
), buf
, output_env_get_xvgr_tlabel(oenv
),
1165 if (output_env_get_print_xvgr_codes(oenv
))
1167 fprintf(fp
, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
1173 if (output_env_get_print_xvgr_codes(oenv
))
1175 fprintf(fp
, "@ subtitle \"after lsq fit to mirror %s\"\n", gn_fit
);
1177 xvgr_legend(fp
, nrms
, (const char**)gn_rms
, oenv
);
1179 for (i
= 0; (i
< teller
); i
++)
1181 if (bSplit
&& i
> 0 && fabs(time
[i
]) < 1e-5)
1183 fprintf(fp
, "%s\n", output_env_get_print_xvgr_codes(oenv
) ? "&" : "");
1185 fprintf(fp
, "%12.7f", time
[i
]);
1186 for (j
= 0; (j
< nrms
); j
++)
1188 fprintf(fp
, " %12.7f", rlsm
[j
][i
]);
1197 sprintf(buf
, "Average %s", whatxvgname
[ewhat
]);
1198 sprintf(buf2
, "Average %s", whatxvglabel
[ewhat
]);
1199 fp
= xvgropen(opt2fn("-a", NFILE
, fnm
), buf
, "Residue", buf2
, oenv
);
1200 for (j
= 0; (j
< nrms
); j
++)
1202 fprintf(fp
, "%10d %10g\n", j
, rlstot
/teller
);
1209 fp
= xvgropen("aver.xvg", gn_rms
[0], "Residue", whatxvglabel
[ewhat
], oenv
);
1210 for (j
= 0; (j
< irms
[0]); j
++)
1212 fprintf(fp
, "%10d %10g\n", j
, rlsnorm
[j
]/teller
);
1216 do_view(oenv
, opt2fn_null("-a", NFILE
, fnm
), "-graphtype bar");
1217 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), NULL
);
1218 do_view(oenv
, opt2fn_null("-mir", NFILE
, fnm
), NULL
);
1219 do_view(oenv
, opt2fn_null("-m", NFILE
, fnm
), NULL
);
1220 do_view(oenv
, opt2fn_null("-bm", NFILE
, fnm
), NULL
);
1221 do_view(oenv
, opt2fn_null("-dist", NFILE
, fnm
), NULL
);