1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #include "gmxpre.h"
39
40 #include <cmath>
41 #include <cstdlib>
42
43 #include <algorithm>
44
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/fileio/confio.h"
48 #include "gromacs/fileio/matio.h"
49 #include "gromacs/fileio/trxio.h"
50 #include "gromacs/fileio/xvgr.h"
51 #include "gromacs/gmxana/cmat.h"
52 #include "gromacs/gmxana/gmx_ana.h"
53 #include "gromacs/gmxana/princ.h"
54 #include "gromacs/math/do_fit.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/topology/ifunc.h"
60 #include "gromacs/topology/index.h"
61 #include "gromacs/topology/topology.h"
62 #include "gromacs/utility/arraysize.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/pleasecite.h"
66 #include "gromacs/utility/smalloc.h"
67
norm_princ(const t_atoms * atoms,int isize,int * index,int natoms,rvec * x)68 static void norm_princ(const t_atoms* atoms, int isize, int* index, int natoms, rvec* x)
69 {
70 int i, m;
71 rvec princ, vec;
72
73 /* equalize principal components: */
74 /* orient principal axes, get principal components */
75 orient_princ(atoms, isize, index, natoms, x, nullptr, princ);
76
77 /* calc our own principal components */
78 clear_rvec(vec);
79 for (m = 0; m < DIM; m++)
80 {
81 for (i = 0; i < isize; i++)
82 {
83 vec[m] += gmx::square(x[index[i]][m]);
84 }
85 vec[m] = std::sqrt(vec[m] / static_cast<real>(isize));
86 /* calculate scaling constants */
87 vec[m] = 1.0 / (std::sqrt(3.0) * vec[m]);
88 }
89
90 /* scale coordinates */
91 for (i = 0; i < natoms; i++)
92 {
93 for (m = 0; m < DIM; m++)
94 {
95 x[i][m] *= vec[m];
96 }
97 }
98 }
99
gmx_rms(int argc,char * argv[])100 int gmx_rms(int argc, char* argv[])
101 {
102 const char* desc[] = {
103 "[THISMODULE] compares two structures by computing the root mean square",
104 "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
105 "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
106 "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
107 "This is selected by [TT]-what[tt].[PAR]",
108
109 "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
110 "reference structure. The reference structure",
111 "is taken from the structure file ([TT]-s[tt]).[PAR]",
112
113 "With option [TT]-mir[tt] also a comparison with the mirror image of",
114 "the reference structure is calculated.",
115 "This is useful as a reference for 'significant' values, see",
116 "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
117
118 "Option [TT]-prev[tt] produces the comparison with a previous frame",
119 "the specified number of frames ago.[PAR]",
120
121 "Option [TT]-m[tt] produces a matrix in [REF].xpm[ref] format of",
122 "comparison values of each structure in the trajectory with respect to",
123 "each other structure. This file can be visualized with for instance",
124 "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]",
125
126 "Option [TT]-fit[tt] controls the least-squares fitting of",
127 "the structures on top of each other: complete fit (rotation and",
128 "translation), translation only, or no fitting at all.[PAR]",
129
130 "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
131 "If you select the option (default) and ",
132 "supply a valid [REF].tpr[ref] file masses will be taken from there, ",
133 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
134 "[TT]GMXLIB[tt]. This is fine for proteins, but not",
135 "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
136 "is assigned to unknown atoms. You can check whether this happened by",
137 "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
138
139 "With [TT]-f2[tt], the 'other structures' are taken from a second",
140 "trajectory, this generates a comparison matrix of one trajectory",
141 "versus the other.[PAR]",
142
143 "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
144
145 "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
146 "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
147 "comparison group are considered."
148 };
149 static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
150 static gmx_bool bDeltaLog = FALSE;
151 static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
152 static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1, bond_user_min = -1,
153 delta_maxy = 0.0;
154 /* strings and things for selecting difference method */
155 enum
156 {
157 ewSel,
158 ewRMSD,
159 ewRho,
160 ewRhoSc,
161 ewNR
162 };
163 int ewhat;
164 const char* what[ewNR + 1] = { nullptr, "rmsd", "rho", "rhosc", nullptr };
165 const char* whatname[ewNR] = { nullptr, "RMSD", "Rho", "Rho sc" };
166 const char* whatlabel[ewNR] = { nullptr, "RMSD (nm)", "Rho", "Rho sc" };
167 const char* whatxvgname[ewNR] = { nullptr, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
168 const char* whatxvglabel[ewNR] = { nullptr, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
169 /* strings and things for fitting methods */
170 enum
171 {
172 efSel,
173 efFit,
174 efReset,
175 efNone,
176 efNR
177 };
178 int efit;
179 const char* fit[efNR + 1] = { nullptr, "rot+trans", "translation", "none", nullptr };
180 const char* fitgraphlabel[efNR + 1] = { nullptr, "lsq fit", "translational fit", "no fit" };
181 static int nrms = 1;
182 static gmx_bool bMassWeighted = TRUE;
183 t_pargs pa[] = {
184 { "-what", FALSE, etENUM, { what }, "Structural difference measure" },
185 { "-pbc", FALSE, etBOOL, { &bPBC }, "PBC check" },
186 { "-fit", FALSE, etENUM, { fit }, "Fit to reference structure" },
187 { "-prev", FALSE, etINT, { &prev }, "Compare with previous frame" },
188 { "-split", FALSE, etBOOL, { &bSplit }, "Split graph where time is zero" },
189 { "-fitall", FALSE, etBOOL, { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
190 { "-skip", FALSE, etINT, { &freq }, "Only write every nr-th frame to matrix" },
191 { "-skip2", FALSE, etINT, { &freq2 }, "Only write every nr-th frame to matrix" },
192 { "-max", FALSE, etREAL, { &rmsd_user_max }, "Maximum level in comparison matrix" },
193 { "-min", FALSE, etREAL, { &rmsd_user_min }, "Minimum level in comparison matrix" },
194 { "-bmax", FALSE, etREAL, { &bond_user_max }, "Maximum level in bond angle matrix" },
195 { "-bmin", FALSE, etREAL, { &bond_user_min }, "Minimum level in bond angle matrix" },
196 { "-mw", FALSE, etBOOL, { &bMassWeighted }, "Use mass weighting for superposition" },
197 { "-nlevels", FALSE, etINT, { &nlevels }, "Number of levels in the matrices" },
198 { "-ng", FALSE, etINT, { &nrms }, "Number of groups to compute RMS between" },
199 { "-dlog", FALSE, etBOOL, { &bDeltaLog }, "HIDDENUse a log x-axis in the delta t matrix" },
200 { "-dmax", FALSE, etREAL, { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
201 { "-aver", FALSE, etINT, { &avl }, "HIDDENAverage over this distance in the RMSD matrix" }
202 };
203 int natoms_trx, natoms_trx2, natoms;
204 int i, j, k, m;
205 #define NFRAME 5000
206 int maxframe = NFRAME, maxframe2 = NFRAME;
207 real t, *w_rls, *w_rms, *w_rls_m = nullptr, *w_rms_m = nullptr;
208 gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
209 gmx_bool bFit, bReset;
210 t_topology top;
211 PbcType pbcType;
212 t_iatom* iatom = nullptr;
213
214 matrix box = { { 0 } };
215 rvec * x, *xp, *xm = nullptr, **mat_x = nullptr, **mat_x2, *mat_x2_j = nullptr, vec1, vec2;
216 t_trxstatus* status;
217 char buf[256], buf2[256];
218 int ncons = 0;
219 FILE* fp;
220 real rlstot = 0, **rls, **rlsm = nullptr, *time, *time2, *rlsnorm = nullptr,
221 **rmsd_mat = nullptr, **bond_mat = nullptr, *axis, *axis2, *del_xaxis, *del_yaxis,
222 rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
223 real ** rmsdav_mat = nullptr, av_tot, weight, weight_tot;
224 real ** delta = nullptr, delta_max, delta_scalex = 0, delta_scaley = 0, *delta_tot;
225 int delta_xsize = 0, del_lev = 100, mx, my, abs_my;
226 gmx_bool bA1, bA2, bPrev, bTop, *bInMat = nullptr;
227 int ifit, *irms, ibond = 0, *ind_bond1 = nullptr, *ind_bond2 = nullptr, n_ind_m = 0;
228 int * ind_fit, **ind_rms, *ind_m = nullptr, *rev_ind_m = nullptr, *ind_rms_m = nullptr;
229 char * gn_fit, **gn_rms;
230 t_rgb rlo, rhi;
231 gmx_output_env_t* oenv;
232 gmx_rmpbc_t gpbc = nullptr;
233
234 t_filenm fnm[] = {
235 { efTPS, nullptr, nullptr, ffREAD }, { efTRX, "-f", nullptr, ffREAD },
236 { efTRX, "-f2", nullptr, ffOPTRD }, { efNDX, nullptr, nullptr, ffOPTRD },
237 { efXVG, nullptr, "rmsd", ffWRITE }, { efXVG, "-mir", "rmsdmir", ffOPTWR },
238 { efXVG, "-a", "avgrp", ffOPTWR }, { efXVG, "-dist", "rmsd-dist", ffOPTWR },
239 { efXPM, "-m", "rmsd", ffOPTWR }, { efDAT, "-bin", "rmsd", ffOPTWR },
240 { efXPM, "-bm", "bond", ffOPTWR }
241 };
242 #define NFILE asize(fnm)
243
244 if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW, NFILE, fnm,
245 asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
246 {
247 return 0;
248 }
249 /* parse enumerated options: */
250 ewhat = nenum(what);
251 if (ewhat == ewRho || ewhat == ewRhoSc)
252 {
253 please_cite(stdout, "Maiorov95");
254 }
255 efit = nenum(fit);
256 bFit = efit == efFit;
257 bReset = efit == efReset;
258 if (bFit)
259 {
260 bReset = TRUE; /* for fit, reset *must* be set */
261 }
262 else
263 {
264 bFitAll = FALSE;
265 }
266
267 /* mark active cmdline options */
268 bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
269 bFile2 = opt2bSet("-f2", NFILE, fnm);
270 bMat = opt2bSet("-m", NFILE, fnm);
271 bBond = opt2bSet("-bm", NFILE, fnm);
272 bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
273 * your RMSD matrix (hidden option */
274 bNorm = opt2bSet("-a", NFILE, fnm);
275 bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
276 if (freq <= 0)
277 {
278 fprintf(stderr,
279 "The number of frames to skip is <= 0. "
280 "Writing out all frames.\n\n");
281 freq = 1;
282 }
283 if (!bFreq2)
284 {
285 freq2 = freq;
286 }
287 else if (bFile2 && freq2 <= 0)
288 {
289 fprintf(stderr,
290 "The number of frames to skip in second trajectory is <= 0.\n"
291 " Writing out all frames.\n\n");
292 freq2 = 1;
293 }
294
295 bPrev = (prev > 0);
296 if (bPrev)
297 {
298 fprintf(stderr,
299 "WARNING: using option -prev with large trajectories will\n"
300 " require a lot of memory and could lead to crashes\n");
301 prev = abs(prev);
302 if (freq != 1)
303 {
304 fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
305 }
306 }
307
308 if (bFile2 && !bMat && !bBond)
309 {
310 fprintf(stderr,
311 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
312 " will not read from %s\n",
313 opt2fn("-f2", NFILE, fnm));
314 bFile2 = FALSE;
315 }
316
317 if (bDelta)
318 {
319 bMat = TRUE;
320 if (bFile2)
321 {
322 fprintf(stderr,
323 "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
324 " will not read from %s\n",
325 opt2fn("-f2", NFILE, fnm));
326 bFile2 = FALSE;
327 }
328 }
329
330 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, TRUE);
331 snew(w_rls, top.atoms.nr);
332 snew(w_rms, top.atoms.nr);
333
334 if (!bTop && bBond)
335 {
336 fprintf(stderr,
337 "WARNING: Need a run input file for bond angle matrix,\n"
338 " will not calculate bond angle matrix.\n");
339 bBond = FALSE;
340 }
341
342 if (bReset)
343 {
344 fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares" : "translational");
345 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit, &ind_fit, &gn_fit);
346 }
347 else
348 {
349 ifit = 0;
350 }
351
352 if (bReset)
353 {
354 if (bFit && ifit < 3)
355 {
356 gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
357 }
358
359 bMass = FALSE;
360 for (i = 0; i < ifit; i++)
361 {
362 if (bMassWeighted)
363 {
364 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
365 }
366 else
367 {
368 w_rls[ind_fit[i]] = 1;
369 }
370 bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
371 }
372 if (!bMass)
373 {
374 fprintf(stderr, "All masses in the fit group are 0, using masses of 1\n");
375 for (i = 0; i < ifit; i++)
376 {
377 w_rls[ind_fit[i]] = 1;
378 }
379 }
380 }
381
382 if (bMat || bBond)
383 {
384 nrms = 1;
385 }
386
387 snew(gn_rms, nrms);
388 snew(ind_rms, nrms);
389 snew(irms, nrms);
390
391 fprintf(stderr, "Select group%s for %s calculation\n", (nrms > 1) ? "s" : "", whatname[ewhat]);
392 get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), nrms, irms, ind_rms, gn_rms);
393
394 if (bNorm)
395 {
396 snew(rlsnorm, irms[0]);
397 }
398 snew(rls, nrms);
399 for (j = 0; j < nrms; j++)
400 {
401 snew(rls[j], maxframe);
402 }
403 if (bMirror)
404 {
405 snew(rlsm, nrms);
406 for (j = 0; j < nrms; j++)
407 {
408 snew(rlsm[j], maxframe);
409 }
410 }
411 snew(time, maxframe);
412 for (j = 0; j < nrms; j++)
413 {
414 bMass = FALSE;
415 for (i = 0; i < irms[j]; i++)
416 {
417 if (bMassWeighted)
418 {
419 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
420 }
421 else
422 {
423 w_rms[ind_rms[j][i]] = 1.0;
424 }
425 bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
426 }
427 if (!bMass)
428 {
429 fprintf(stderr, "All masses in group %d are 0, using masses of 1\n", j);
430 for (i = 0; i < irms[j]; i++)
431 {
432 w_rms[ind_rms[j][i]] = 1;
433 }
434 }
435 }
436 /* Prepare reference frame */
437 if (bPBC)
438 {
439 gpbc = gmx_rmpbc_init(&top.idef, pbcType, top.atoms.nr);
440 gmx_rmpbc(gpbc, top.atoms.nr, box, xp);
441 }
442 if (bReset)
443 {
444 reset_x(ifit, ind_fit, top.atoms.nr, nullptr, xp, w_rls);
445 }
446 if (bMirror)
447 {
448 /* generate reference structure mirror image: */
449 snew(xm, top.atoms.nr);
450 for (i = 0; i < top.atoms.nr; i++)
451 {
452 copy_rvec(xp[i], xm[i]);
453 xm[i][XX] = -xm[i][XX];
454 }
455 }
456 if (ewhat == ewRhoSc)
457 {
458 norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
459 }
460
461 /* read first frame */
462 natoms_trx = read_first_x(oenv, &status, opt2fn("-f", NFILE, fnm), &t, &x, box);
463 if (natoms_trx != top.atoms.nr)
464 {
465 fprintf(stderr, "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
466 top.atoms.nr, natoms_trx);
467 }
468 natoms = std::min(top.atoms.nr, natoms_trx);
469 if (bMat || bBond || bPrev)
470 {
471 snew(mat_x, NFRAME);
472
473 if (bPrev)
474 {
475 /* With -prev we use all atoms for simplicity */
476 n_ind_m = natoms;
477 }
478 else
479 {
480 /* Check which atoms we need (fit/rms) */
481 snew(bInMat, natoms);
482 for (i = 0; i < ifit; i++)
483 {
484 bInMat[ind_fit[i]] = TRUE;
485 }
486 n_ind_m = ifit;
487 for (i = 0; i < irms[0]; i++)
488 {
489 if (!bInMat[ind_rms[0][i]])
490 {
491 bInMat[ind_rms[0][i]] = TRUE;
492 n_ind_m++;
493 }
494 }
495 }
496 /* Make an index of needed atoms */
497 snew(ind_m, n_ind_m);
498 snew(rev_ind_m, natoms);
499 j = 0;
500 for (i = 0; i < natoms; i++)
501 {
502 if (bPrev || bInMat[i])
503 {
504 ind_m[j] = i;
505 rev_ind_m[i] = j;
506 j++;
507 }
508 }
509 snew(w_rls_m, n_ind_m);
510 snew(ind_rms_m, irms[0]);
511 snew(w_rms_m, n_ind_m);
512 for (i = 0; i < ifit; i++)
513 {
514 w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
515 }
516 for (i = 0; i < irms[0]; i++)
517 {
518 ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
519 w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
520 }
521 sfree(bInMat);
522 }
523
524 if (bBond)
525 {
526 ncons = 0;
527 for (k = 0; k < F_NRE; k++)
528 {
529 if (IS_CHEMBOND(k))
530 {
531 ncons += top.idef.il[k].nr / 3;
532 }
533 }
534 fprintf(stderr, "Found %d bonds in topology\n", ncons);
535 snew(ind_bond1, ncons);
536 snew(ind_bond2, ncons);
537 ibond = 0;
538 for (k = 0; k < F_NRE; k++)
539 {
540 if (IS_CHEMBOND(k))
541 {
542 iatom = top.idef.il[k].iatoms;
543 ncons = top.idef.il[k].nr / 3;
544 for (i = 0; i < ncons; i++)
545 {
546 bA1 = FALSE;
547 bA2 = FALSE;
548 for (j = 0; j < irms[0]; j++)
549 {
550 if (iatom[3 * i + 1] == ind_rms[0][j])
551 {
552 bA1 = TRUE;
553 }
554 if (iatom[3 * i + 2] == ind_rms[0][j])
555 {
556 bA2 = TRUE;
557 }
558 }
559 if (bA1 && bA2)
560 {
561 ind_bond1[ibond] = rev_ind_m[iatom[3 * i + 1]];
562 ind_bond2[ibond] = rev_ind_m[iatom[3 * i + 2]];
563 ibond++;
564 }
565 }
566 }
567 }
568 fprintf(stderr, "Using %d bonds for bond angle matrix\n", ibond);
569 if (ibond == 0)
570 {
571 gmx_fatal(FARGS, "0 bonds found");
572 }
573 }
574
575 /* start looping over frames: */
576 int tel_mat = 0;
577 int teller = 0;
578 int frame = 0;
579 do
580 {
581 if (bPBC)
582 {
583 gmx_rmpbc(gpbc, natoms, box, x);
584 }
585
586 if (bReset)
587 {
588 reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
589 }
590 if (ewhat == ewRhoSc)
591 {
592 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
593 }
594
595 if (bFit)
596 {
597 /*do the least squares fit to original structure*/
598 do_fit(natoms, w_rls, xp, x);
599 }
600
601 if (frame % freq == 0)
602 {
603 /* keep frame for matrix calculation */
604 if (bMat || bBond || bPrev)
605 {
606 if (tel_mat >= NFRAME)
607 {
608 srenew(mat_x, tel_mat + 1);
609 }
610 snew(mat_x[tel_mat], n_ind_m);
611 for (i = 0; i < n_ind_m; i++)
612 {
613 copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]);
614 }
615 }
616 tel_mat++;
617
618 /*calculate energy of root_least_squares*/
619 if (bPrev)
620 {
621 j = tel_mat - prev - 1;
622 if (j < 0)
623 {
624 j = 0;
625 }
626 for (i = 0; i < n_ind_m; i++)
627 {
628 copy_rvec(mat_x[j][i], xp[ind_m[i]]);
629 }
630 if (bReset)
631 {
632 reset_x(ifit, ind_fit, natoms, nullptr, xp, w_rls);
633 }
634 if (bFit)
635 {
636 do_fit(natoms, w_rls, x, xp);
637 }
638 }
639 for (j = 0; (j < nrms); j++)
640 {
641 rls[j][teller] = calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
642 }
643 if (bNorm)
644 {
645 for (j = 0; (j < irms[0]); j++)
646 {
647 rlsnorm[j] += calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp);
648 }
649 }
650
651 if (bMirror)
652 {
653 if (bFit)
654 {
655 /*do the least squares fit to mirror of original structure*/
656 do_fit(natoms, w_rls, xm, x);
657 }
658
659 for (j = 0; j < nrms; j++)
660 {
661 rlsm[j][teller] =
662 calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm);
663 }
664 }
665 time[teller] = output_env_conv_time(oenv, t);
666
667 teller++;
668 }
669 frame++;
670 if (teller >= maxframe)
671 {
672 maxframe += NFRAME;
673 srenew(time, maxframe);
674 for (j = 0; (j < nrms); j++)
675 {
676 srenew(rls[j], maxframe);
677 }
678 if (bMirror)
679 {
680 for (j = 0; (j < nrms); j++)
681 {
682 srenew(rlsm[j], maxframe);
683 }
684 }
685 }
686 } while (read_next_x(oenv, status, &t, x, box));
687 close_trx(status);
688
689 int tel_mat2 = 0;
690 int teller2 = 0;
691 int frame2 = 0;
692 if (bFile2)
693 {
694 snew(time2, maxframe2);
695
696 fprintf(stderr, "\nWill read second trajectory file\n");
697 snew(mat_x2, NFRAME);
698 natoms_trx2 = read_first_x(oenv, &status, opt2fn("-f2", NFILE, fnm), &t, &x, box);
699 if (natoms_trx2 != natoms_trx)
700 {
701 gmx_fatal(FARGS,
702 "Second trajectory (%d atoms) does not match the first one"
703 " (%d atoms)",
704 natoms_trx2, natoms_trx);
705 }
706 frame2 = 0;
707 do
708 {
709 if (bPBC)
710 {
711 gmx_rmpbc(gpbc, natoms, box, x);
712 }
713
714 if (bReset)
715 {
716 reset_x(ifit, ind_fit, natoms, nullptr, x, w_rls);
717 }
718 if (ewhat == ewRhoSc)
719 {
720 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
721 }
722
723 if (bFit)
724 {
725 /*do the least squares fit to original structure*/
726 do_fit(natoms, w_rls, xp, x);
727 }
728
729 if (frame2 % freq2 == 0)
730 {
731 /* keep frame for matrix calculation */
732 if (bMat)
733 {
734 if (tel_mat2 >= NFRAME)
735 {
736 srenew(mat_x2, tel_mat2 + 1);
737 }
738 snew(mat_x2[tel_mat2], n_ind_m);
739 for (i = 0; i < n_ind_m; i++)
740 {
741 copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]);
742 }
743 }
744 tel_mat2++;
745
746 time2[teller2] = output_env_conv_time(oenv, t);
747
748 teller2++;
749 }
750 frame2++;
751 if (teller2 >= maxframe2)
752 {
753 maxframe2 += NFRAME;
754 srenew(time2, maxframe2);
755 }
756 } while (read_next_x(oenv, status, &t, x, box));
757 close_trx(status);
758 }
759 else
760 {
761 mat_x2 = mat_x;
762 time2 = time;
763 tel_mat2 = tel_mat;
764 freq2 = freq;
765 }
766 gmx_rmpbc_done(gpbc);
767
768 if (bMat || bBond)
769 {
770 /* calculate RMS matrix */
771 fprintf(stderr, "\n");
772 if (bMat)
773 {
774 fprintf(stderr, "Building %s matrix, %dx%d elements\n", whatname[ewhat], tel_mat, tel_mat2);
775 snew(rmsd_mat, tel_mat);
776 }
777 if (bBond)
778 {
779 fprintf(stderr, "Building bond angle matrix, %dx%d elements\n", tel_mat, tel_mat2);
780 snew(bond_mat, tel_mat);
781 }
782 snew(axis, tel_mat);
783 snew(axis2, tel_mat2);
784 rmsd_max = 0;
785 if (bFile2)
786 {
787 rmsd_min = 1e10;
788 }
789 else
790 {
791 rmsd_min = 0;
792 }
793 rmsd_avg = 0;
794 bond_max = 0;
795 bond_min = 1e10;
796 for (j = 0; j < tel_mat2; j++)
797 {
798 axis2[j] = time2[freq2 * j];
799 }
800 if (bDelta)
801 {
802 if (bDeltaLog)
803 {
804 delta_scalex = 8.0 / std::log(2.0);
805 delta_xsize = gmx::roundToInt(std::log(tel_mat / 2.) * delta_scalex) + 1;
806 }
807 else
808 {
809 delta_xsize = tel_mat / 2;
810 }
811 delta_scaley = 1.0 / delta_maxy;
812 snew(delta, delta_xsize);
813 for (j = 0; j < delta_xsize; j++)
814 {
815 snew(delta[j], del_lev + 1);
816 }
817 if (avl > 0)
818 {
819 snew(rmsdav_mat, tel_mat);
820 for (j = 0; j < tel_mat; j++)
821 {
822 snew(rmsdav_mat[j], tel_mat);
823 }
824 }
825 }
826
827 if (bFitAll)
828 {
829 snew(mat_x2_j, natoms);
830 }
831 for (i = 0; i < tel_mat; i++)
832 {
833 axis[i] = time[freq * i];
834 fprintf(stderr, "\r element %5d; time %5.2f ", i, axis[i]);
835 fflush(stderr);
836 if (bMat)
837 {
838 snew(rmsd_mat[i], tel_mat2);
839 }
840 if (bBond)
841 {
842 snew(bond_mat[i], tel_mat2);
843 }
844 for (j = 0; j < tel_mat2; j++)
845 {
846 if (bFitAll)
847 {
848 for (k = 0; k < n_ind_m; k++)
849 {
850 copy_rvec(mat_x2[j][k], mat_x2_j[k]);
851 }
852 do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j);
853 }
854 else
855 {
856 mat_x2_j = mat_x2[j];
857 }
858 if (bMat)
859 {
860 if (bFile2 || (i < j))
861 {
862 rmsd_mat[i][j] = calc_similar_ind(ewhat != ewRMSD, irms[0], ind_rms_m,
863 w_rms_m, mat_x[i], mat_x2_j);
864 if (rmsd_mat[i][j] > rmsd_max)
865 {
866 rmsd_max = rmsd_mat[i][j];
867 }
868 if (rmsd_mat[i][j] < rmsd_min)
869 {
870 rmsd_min = rmsd_mat[i][j];
871 }
872 rmsd_avg += rmsd_mat[i][j];
873 }
874 else
875 {
876 rmsd_mat[i][j] = rmsd_mat[j][i];
877 }
878 }
879 if (bBond)
880 {
881 if (bFile2 || (i <= j))
882 {
883 ang = 0.0;
884 for (m = 0; m < ibond; m++)
885 {
886 rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1);
887 rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2);
888 ang += std::acos(cos_angle(vec1, vec2));
889 }
890 bond_mat[i][j] = ang * 180.0 / (M_PI * ibond);
891 if (bond_mat[i][j] > bond_max)
892 {
893 bond_max = bond_mat[i][j];
894 }
895 if (bond_mat[i][j] < bond_min)
896 {
897 bond_min = bond_mat[i][j];
898 }
899 }
900 else
901 {
902 bond_mat[i][j] = bond_mat[j][i];
903 }
904 }
905 }
906 }
907 if (bFile2)
908 {
909 rmsd_avg /= static_cast<real>(tel_mat) * static_cast<real>(tel_mat2);
910 }
911 else
912 {
913 rmsd_avg /= tel_mat * (tel_mat - 1) / 2.;
914 }
915 if (bMat && (avl > 0))
916 {
917 rmsd_max = 0.0;
918 rmsd_min = 0.0;
919 rmsd_avg = 0.0;
920 for (j = 0; j < tel_mat - 1; j++)
921 {
922 for (i = j + 1; i < tel_mat; i++)
923 {
924 av_tot = 0;
925 weight_tot = 0;
926 for (my = -avl; my <= avl; my++)
927 {
928 if ((j + my >= 0) && (j + my < tel_mat))
929 {
930 abs_my = std::abs(my);
931 for (mx = -avl; mx <= avl; mx++)
932 {
933 if ((i + mx >= 0) && (i + mx < tel_mat))
934 {
935 weight = avl + 1.0 - std::max(std::abs(mx), abs_my);
936 av_tot += weight * rmsd_mat[i + mx][j + my];
937 weight_tot += weight;
938 }
939 }
940 }
941 }
942 rmsdav_mat[i][j] = av_tot / weight_tot;
943 rmsdav_mat[j][i] = rmsdav_mat[i][j];
944 if (rmsdav_mat[i][j] > rmsd_max)
945 {
946 rmsd_max = rmsdav_mat[i][j];
947 }
948 }
949 }
950 rmsd_mat = rmsdav_mat;
951 }
952
953 if (bMat)
954 {
955 fprintf(stderr, "\n%s: Min %f, Max %f, Avg %f\n", whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg);
956 rlo.r = 1;
957 rlo.g = 1;
958 rlo.b = 1;
959 rhi.r = 0;
960 rhi.g = 0;
961 rhi.b = 0;
962 if (rmsd_user_max != -1)
963 {
964 rmsd_max = rmsd_user_max;
965 }
966 if (rmsd_user_min != -1)
967 {
968 rmsd_min = rmsd_user_min;
969 }
970 if ((rmsd_user_max != -1) || (rmsd_user_min != -1))
971 {
972 fprintf(stderr, "Min and Max value set to resp. %f and %f\n", rmsd_min, rmsd_max);
973 }
974 sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]);
975 write_xpm(opt2FILE("-m", NFILE, fnm, "w"), 0, buf, whatlabel[ewhat],
976 output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat,
977 tel_mat2, axis, axis2, rmsd_mat, rmsd_min, rmsd_max, rlo, rhi, &nlevels);
978 /* Print the distribution of RMSD values */
979 if (opt2bSet("-dist", NFILE, fnm))
980 {
981 low_rmsd_dist(opt2fn("-dist", NFILE, fnm), rmsd_max, tel_mat, rmsd_mat, oenv);
982 }
983
984 if (bDelta)
985 {
986 snew(delta_tot, delta_xsize);
987 for (j = 0; j < tel_mat - 1; j++)
988 {
989 for (i = j + 1; i < tel_mat; i++)
990 {
991 mx = i - j;
992 if (mx < tel_mat / 2)
993 {
994 if (bDeltaLog)
995 {
996 mx = gmx::roundToInt(std::log(static_cast<real>(mx)) * delta_scalex);
997 }
998 my = gmx::roundToInt(rmsd_mat[i][j] * delta_scaley * static_cast<real>(del_lev));
999 delta_tot[mx] += 1.0;
1000 if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy))
1001 {
1002 delta[mx][my] += 1.0;
1003 }
1004 }
1005 }
1006 }
1007 delta_max = 0;
1008 for (i = 0; i < delta_xsize; i++)
1009 {
1010 if (delta_tot[i] > 0.0)
1011 {
1012 delta_tot[i] = 1.0 / delta_tot[i];
1013 for (j = 0; j <= del_lev; j++)
1014 {
1015 delta[i][j] *= delta_tot[i];
1016 if (delta[i][j] > delta_max)
1017 {
1018 delta_max = delta[i][j];
1019 }
1020 }
1021 }
1022 }
1023 fprintf(stderr, "Maximum in delta matrix: %f\n", delta_max);
1024 snew(del_xaxis, delta_xsize);
1025 snew(del_yaxis, del_lev + 1);
1026 for (i = 0; i < delta_xsize; i++)
1027 {
1028 del_xaxis[i] = axis[i] - axis[0];
1029 }
1030 for (i = 0; i < del_lev + 1; i++)
1031 {
1032 del_yaxis[i] = delta_maxy * static_cast<real>(i) / static_cast<real>(del_lev);
1033 }
1034 sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]);
1035 fp = gmx_ffopen("delta.xpm", "w");
1036 write_xpm(fp, 0, buf, "density", output_env_get_time_label(oenv), whatlabel[ewhat],
1037 delta_xsize, del_lev + 1, del_xaxis, del_yaxis, delta, 0.0, delta_max,
1038 rlo, rhi, &nlevels);
1039 gmx_ffclose(fp);
1040 }
1041 if (opt2bSet("-bin", NFILE, fnm))
1042 {
1043 /* NB: File must be binary if we use fwrite */
1044 fp = ftp2FILE(efDAT, NFILE, fnm, "wb");
1045 for (i = 0; i < tel_mat; i++)
1046 {
1047 if (static_cast<int>(fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp)) != tel_mat2)
1048 {
1049 gmx_fatal(FARGS, "Error writing to output file");
1050 }
1051 }
1052 gmx_ffclose(fp);
1053 }
1054 }
1055 if (bBond)
1056 {
1057 fprintf(stderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max);
1058 if (bond_user_max != -1)
1059 {
1060 bond_max = bond_user_max;
1061 }
1062 if (bond_user_min != -1)
1063 {
1064 bond_min = bond_user_min;
1065 }
1066 if ((bond_user_max != -1) || (bond_user_min != -1))
1067 {
1068 fprintf(stderr,
1069 "Bond angle Min and Max set to:\n"
1070 "Min. angle: %f, Max. angle: %f\n",
1071 bond_min, bond_max);
1072 }
1073 rlo.r = 1;
1074 rlo.g = 1;
1075 rlo.b = 1;
1076 rhi.r = 0;
1077 rhi.g = 0;
1078 rhi.b = 0;
1079 sprintf(buf, "%s av. bond angle deviation", gn_rms[0]);
1080 write_xpm(opt2FILE("-bm", NFILE, fnm, "w"), 0, buf, "degrees",
1081 output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat,
1082 tel_mat2, axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels);
1083 }
1084 }
1085
1086 bAv = opt2bSet("-a", NFILE, fnm);
1087
1088 /* Write the RMSD's to file */
1089 if (!bPrev)
1090 {
1091 sprintf(buf, "%s", whatxvgname[ewhat]);
1092 }
1093 else
1094 {
1095 sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat], time[prev * freq] - time[0],
1096 output_env_get_time_label(oenv).c_str());
1097 }
1098 fp = xvgropen(opt2fn("-o", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv),
1099 whatxvglabel[ewhat], oenv);
1100 if (output_env_get_print_xvgr_codes(oenv))
1101 {
1102 fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n", (nrms == 1) ? "" : "of ", gn_rms[0],
1103 fitgraphlabel[efit], bFit ? " to " : "", bFit ? gn_fit : "");
1104 }
1105 if (nrms != 1)
1106 {
1107 xvgr_legend(fp, nrms, gn_rms, oenv);
1108 }
1109 for (i = 0; (i < teller); i++)
1110 {
1111 if (bSplit && i > 0 && std::abs(time[bPrev ? freq * i : i] / output_env_get_time_factor(oenv)) < 1e-5)
1112 {
1113 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1114 }
1115 fprintf(fp, "%12.7f", time[bPrev ? freq * i : i]);
1116 for (j = 0; (j < nrms); j++)
1117 {
1118 fprintf(fp, " %12.7f", rls[j][i]);
1119 if (bAv)
1120 {
1121 rlstot += rls[j][i];
1122 }
1123 }
1124 fprintf(fp, "\n");
1125 }
1126 xvgrclose(fp);
1127
1128 if (bMirror)
1129 {
1130 /* Write the mirror RMSD's to file */
1131 sprintf(buf, "%s with Mirror", whatxvgname[ewhat]);
1132 sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]);
1133 fp = xvgropen(opt2fn("-mir", NFILE, fnm), buf, output_env_get_xvgr_tlabel(oenv), buf2, oenv);
1134 if (nrms == 1)
1135 {
1136 if (output_env_get_print_xvgr_codes(oenv))
1137 {
1138 fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n", gn_rms[0],
1139 bFit ? gn_fit : "");
1140 }
1141 }
1142 else
1143 {
1144 if (output_env_get_print_xvgr_codes(oenv))
1145 {
1146 fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", bFit ? gn_fit : "");
1147 }
1148 xvgr_legend(fp, nrms, gn_rms, oenv);
1149 }
1150 for (i = 0; (i < teller); i++)
1151 {
1152 if (bSplit && i > 0 && std::abs(time[i]) < 1e-5)
1153 {
1154 fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
1155 }
1156 fprintf(fp, "%12.7f", time[i]);
1157 for (j = 0; (j < nrms); j++)
1158 {
1159 fprintf(fp, " %12.7f", rlsm[j][i]);
1160 }
1161 fprintf(fp, "\n");
1162 }
1163 xvgrclose(fp);
1164 }
1165
1166 if (bAv)
1167 {
1168 sprintf(buf, "Average %s", whatxvgname[ewhat]);
1169 sprintf(buf2, "Average %s", whatxvglabel[ewhat]);
1170 fp = xvgropen(opt2fn("-a", NFILE, fnm), buf, "Residue", buf2, oenv);
1171 for (j = 0; (j < nrms); j++)
1172 {
1173 fprintf(fp, "%10d %10g\n", j, rlstot / static_cast<real>(teller));
1174 }
1175 xvgrclose(fp);
1176 }
1177
1178 if (bNorm)
1179 {
1180 fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv);
1181 for (j = 0; (j < irms[0]); j++)
1182 {
1183 fprintf(fp, "%10d %10g\n", j, rlsnorm[j] / static_cast<real>(teller));
1184 }
1185 xvgrclose(fp);
1186 }
1187 do_view(oenv, opt2fn_null("-a", NFILE, fnm), "-graphtype bar");
1188 do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr);
1189 do_view(oenv, opt2fn_null("-mir", NFILE, fnm), nullptr);
1190 do_view(oenv, opt2fn_null("-m", NFILE, fnm), nullptr);
1191 do_view(oenv, opt2fn_null("-bm", NFILE, fnm), nullptr);
1192 do_view(oenv, opt2fn_null("-dist", NFILE, fnm), nullptr);
1193
1194 return 0;
1195 }
1196