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,2017,2018,2019, 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.
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/commandline/viewit.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/matio.h"
48 #include "gromacs/fileio/trxio.h"
49 #include "gromacs/fileio/xvgr.h"
50 #include "gromacs/gmxana/gmx_ana.h"
51 #include "gromacs/math/functions.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/pbcutil/pbc.h"
55 #include "gromacs/topology/index.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/utility/arraysize.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
62 #include "gromacs/utility/strdb.h"
65 static void calc_dist(int nind, const int index[], const rvec x[], int ePBC, matrix box, real** d)
72 set_pbc(&pbc, ePBC, box);
73 for (i = 0; (i < nind - 1); i++)
75 const real* xi = x[index[i]];
76 for (j = i + 1; (j < nind); j++)
78 pbc_dx(&pbc, xi, x[index[j]], dx);
80 d[i][j] = std::sqrt(temp2);
85 static void calc_dist_tot(int nind,
99 real temp, temp2, temp1_3;
103 set_pbc(&pbc, ePBC, box);
104 for (i = 0; (i < nind - 1); i++)
107 for (j = i + 1; (j < nind); j++)
109 pbc_dx(&pbc, xi, x[index[j]], dx);
111 temp = std::sqrt(temp2);
114 dtot2[i][j] += temp2;
117 temp1_3 = 1.0 / (temp * temp2);
118 dtot1_3[i][j] += temp1_3;
119 dtot1_6[i][j] += temp1_3 * temp1_3;
125 static void calc_nmr(int nind, int nframes, real** dtot1_3, real** dtot1_6, real* max1_3, real* max1_6)
128 real temp1_3, temp1_6;
130 for (i = 0; (i < nind - 1); i++)
132 for (j = i + 1; (j < nind); j++)
134 temp1_3 = gmx::invcbrt(dtot1_3[i][j] / static_cast<real>(nframes));
135 temp1_6 = gmx::invsixthroot(dtot1_6[i][j] / static_cast<real>(nframes));
136 if (temp1_3 > *max1_3)
140 if (temp1_6 > *max1_6)
144 dtot1_3[i][j] = temp1_3;
145 dtot1_6[i][j] = temp1_6;
146 dtot1_3[j][i] = temp1_3;
147 dtot1_6[j][i] = temp1_6;
152 static char Hnum[] = "123";
181 static int read_equiv(const char* eq_fn, t_equiv*** equivptr)
184 char line[STRLEN], resname[10], atomname[10], *lp;
185 int neq, na, n, resnr;
188 fp = gmx_ffopen(eq_fn, "r");
191 while (get_a_line(fp, line, STRLEN))
194 /* this is not efficient, but I'm lazy */
195 srenew(equiv, neq + 1);
196 equiv[neq] = nullptr;
198 if (sscanf(lp, "%s %n", atomname, &n) == 1)
202 equiv[neq][0].nname = gmx_strdup(atomname);
203 while (sscanf(lp, "%d %s %s %n", &resnr, resname, atomname, &n) == 3)
205 /* this is not efficient, but I'm lazy (again) */
206 srenew(equiv[neq], na + 1);
207 equiv[neq][na].set = true;
208 equiv[neq][na].rnr = resnr - 1;
209 equiv[neq][na].rname = gmx_strdup(resname);
210 equiv[neq][na].aname = gmx_strdup(atomname);
213 equiv[neq][na].nname = nullptr;
219 /* make empty element as flag for end of array */
220 srenew(equiv[neq], na + 1);
221 equiv[neq][na].set = false;
222 equiv[neq][na].rnr = 0;
223 equiv[neq][na].rname = nullptr;
224 equiv[neq][na].aname = nullptr;
236 static void dump_equiv(FILE* out, int neq, t_equiv** equiv)
240 fprintf(out, "Dumping equivalent list\n");
241 for (i = 0; i < neq; i++)
243 fprintf(out, "%s", equiv[i][0].nname);
244 for (j = 0; equiv[i][j].set; j++)
246 fprintf(out, " %d %s %s", equiv[i][j].rnr, equiv[i][j].rname, equiv[i][j].aname);
253 is_equiv(int neq, t_equiv** equiv, char** nname, int rnr1, char* rname1, char* aname1, int rnr2, char* rname2, char* aname2)
259 /* we can terminate each loop when bFound is true! */
260 for (i = 0; i < neq && !bFound; i++)
262 /* find first atom */
263 for (j = 0; equiv[i][j].set && !bFound; j++)
265 bFound = (equiv[i][j].rnr == rnr1 && std::strcmp(equiv[i][j].rname, rname1) == 0
266 && std::strcmp(equiv[i][j].aname, aname1) == 0);
270 /* find second atom */
272 for (j = 0; equiv[i][j].set && !bFound; j++)
274 bFound = (equiv[i][j].rnr == rnr2 && std::strcmp(equiv[i][j].rname, rname2) == 0
275 && std::strcmp(equiv[i][j].aname, aname2) == 0);
281 *nname = gmx_strdup(equiv[i - 1][0].nname);
287 static int analyze_noe_equivalent(const char* eq_fn,
288 const t_atoms* atoms,
295 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
296 char * anmi, *anmj, **nnm;
297 gmx_bool bMatch, bEquiv;
305 neq = read_equiv(eq_fn, &equiv);
308 dump_equiv(debug, neq, equiv);
318 for (i = 0; i < isize; i++)
320 if (equiv && i < isize - 1)
322 /* check explicit list of equivalent atoms */
326 rnri = atoms->atom[index[i]].resind;
327 rnrj = atoms->atom[index[j]].resind;
328 bEquiv = is_equiv(neq, equiv, &nnm[i], rnri, *atoms->resinfo[rnri].name,
329 *atoms->atomname[index[i]], rnrj, *atoms->resinfo[rnrj].name,
330 *atoms->atomname[index[j]]);
331 if (nnm[i] && bEquiv)
333 nnm[j] = gmx_strdup(nnm[i]);
337 /* set index for matching atom */
338 noe_index[i] = groupnr;
339 /* skip matching atom */
342 } while (bEquiv && i < isize - 1);
350 /* look for triplets of consecutive atoms with name XX?,
351 X are any number of letters or digits and ? goes from 1 to 3
352 This is supposed to cover all CH3 groups and the like */
353 anmi = *atoms->atomname[index[i]];
354 anmil = std::strlen(anmi);
355 bMatch = i <= isize - 3 && anmi[anmil - 1] == '1';
358 for (j = 1; j < 3; j++)
360 anmj = *atoms->atomname[index[i + j]];
361 anmjl = std::strlen(anmj);
363 && (anmil == anmjl && anmj[anmjl - 1] == Hnum[j]
364 && std::strncmp(anmi, anmj, anmil - 1) == 0);
367 /* set index for this atom */
368 noe_index[i] = groupnr;
371 /* set index for next two matching atoms */
372 for (j = 1; j < 3; j++)
374 noe_index[i + j] = groupnr;
376 /* skip matching atoms */
385 /* make index without looking for equivalent atoms */
386 for (i = 0; i < isize; i++)
392 noe_index[isize] = groupnr;
397 for (i = 0; i < isize; i++)
399 rnri = atoms->atom[index[i]].resind;
400 fprintf(debug, "%s %s %d -> %s\n", *atoms->atomname[index[i]],
401 *atoms->resinfo[rnri].name, rnri, nnm[i] ? nnm[i] : "");
405 for (i = 0; i < isize; i++)
408 if (!noe_gr[gi].aname)
411 noe_gr[gi].anr = index[i];
414 noe_gr[gi].aname = gmx_strdup(nnm[i]);
418 noe_gr[gi].aname = gmx_strdup(*atoms->atomname[index[i]]);
419 if (noe_index[i] == noe_index[i + 1])
421 noe_gr[gi].aname[std::strlen(noe_gr[gi].aname) - 1] = '*';
424 noe_gr[gi].rnr = atoms->atom[index[i]].resind;
425 noe_gr[gi].rname = gmx_strdup(*atoms->resinfo[noe_gr[gi].rnr].name);
426 /* dump group definitions */
429 fprintf(debug, "%d %d %d %d %s %s %d\n", i, gi, noe_gr[gi].ianr, noe_gr[gi].anr,
430 noe_gr[gi].aname, noe_gr[gi].rname, noe_gr[gi].rnr);
434 for (i = 0; i < isize; i++)
443 /* #define NSCALE 3 */
444 /* static char *noe_scale[NSCALE+1] = { */
445 /* "strong", "medium", "weak", "none" */
449 static char* noe2scale(real r3, real r6, real rmax)
452 static char buf[NSCALE + 1];
454 /* r goes from 0 to rmax
455 NSCALE*r/rmax goes from 0 to NSCALE
456 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
457 s3 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE * r3 / rmax));
458 s6 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE * r6 / rmax));
460 for (i = 0; i < s3; i++)
473 static void calc_noe(int isize, const int* noe_index, real** dtot1_3, real** dtot1_6, int gnr, t_noe** noe)
477 /* make half matrix of noe-group distances from atom distances */
478 for (i = 0; i < isize; i++)
481 for (j = i; j < isize; j++)
485 noe[gi][gj].i_3 += 1.0 / gmx::power3(dtot1_3[i][j]);
486 noe[gi][gj].i_6 += 1.0 / gmx::power6(dtot1_6[i][j]);
491 for (i = 0; i < gnr; i++)
493 for (j = i + 1; j < gnr; j++)
495 noe[i][j].r_3 = gmx::invcbrt(noe[i][j].i_3 / static_cast<real>(noe[i][j].nr));
496 noe[i][j].r_6 = gmx::invsixthroot(noe[i][j].i_6 / static_cast<real>(noe[i][j].nr));
497 noe[j][i] = noe[i][j];
502 static void write_noe(FILE* fp, int gnr, t_noe** noe, t_noe_gr* noe_gr, real rmax)
505 real r3, r6, min3, min6;
506 char buf[10], b3[10], b6[10];
510 fprintf(fp, ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n", "ianr", "anr",
511 "anm", "rnm", "rnr", "ianr", "anr", "anm", "rnm", "rnr", "1/r^3", "1/r^6", "intnsty",
512 "Dr", "Da", "scale");
513 for (i = 0; i < gnr; i++)
516 for (j = i + 1; j < gnr; j++)
521 min3 = std::min(r3, min3);
522 min6 = std::min(r6, min6);
523 if (r3 < rmax || r6 < rmax)
525 if (grj.rnr == gri.rnr)
527 sprintf(buf, "%2d", grj.anr - gri.anr);
535 sprintf(b3, "%-5.3f", r3);
539 std::strcpy(b3, "-");
543 sprintf(b6, "%-5.3f", r6);
547 std::strcpy(b6, "-");
549 fprintf(fp, "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
550 gri.ianr + 1, gri.anr + 1, gri.aname, gri.rname, gri.rnr + 1, grj.ianr + 1,
551 grj.anr + 1, grj.aname, grj.rname, grj.rnr + 1, b3, b6,
552 gmx::roundToInt(noe[i][j].i_6), grj.rnr - gri.rnr, buf, noe2scale(r3, r6, rmax));
556 #define MINI ((i == 3) ? min3 : min6)
557 for (i = 3; i <= 6; i += 3)
562 "NOTE: no 1/r^%d averaged distances found below %g, "
568 fprintf(stdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI);
574 static void calc_rms(int nind,
586 real mean, mean2, rms, rmsc;
587 /* N.B. dtot and dtot2 contain the total distance and the total squared
588 * distance respectively, BUT they return RMS and the scaled RMS resp.
595 for (i = 0; (i < nind - 1); i++)
597 for (j = i + 1; (j < nind); j++)
599 mean = dtot[i][j] / static_cast<real>(nframes);
600 mean2 = dtot2[i][j] / static_cast<real>(nframes);
601 rms = std::sqrt(std::max(0.0_real, mean2 - mean * mean));
615 meanmat[i][j] = meanmat[j][i] = mean;
616 rmsmat[i][j] = rmsmat[j][i] = rms;
617 rmscmat[i][j] = rmscmat[j][i] = rmsc;
622 static real rms_diff(int natom, real** d, real** d_r)
628 for (i = 0; (i < natom - 1); i++)
630 for (j = i + 1; (j < natom); j++)
632 r = d[i][j] - d_r[i][j];
636 r2 /= gmx::exactDiv(natom * (natom - 1), 2);
638 return std::sqrt(r2);
641 int gmx_rmsdist(int argc, char* argv[])
643 const char* desc[] = {
644 "[THISMODULE] computes the root mean square deviation of atom distances,",
645 "which has the advantage that no fit is needed like in standard RMS",
646 "deviation as computed by [gmx-rms].",
647 "The reference structure is taken from the structure file.",
648 "The RMSD at time t is calculated as the RMS",
649 "of the differences in distance between atom-pairs in the reference",
650 "structure and the structure at time t.[PAR]",
651 "[THISMODULE] can also produce matrices of the rms distances, rms distances",
652 "scaled with the mean distance and the mean distances and matrices with",
653 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
654 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
655 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
656 "can be generated, by default averaging over equivalent hydrogens",
657 "(all triplets of hydrogens named \\*[123]). Additionally a list of",
658 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
659 "a set of equivalent atoms specified as residue number and name and",
660 "atom name; e.g.:[PAR]",
661 "[TT]HB* 3 SER HB1 3 SER HB2[tt][PAR]",
662 "Residue and atom names must exactly match those in the structure",
663 "file, including case. Specifying non-sequential atoms is undefined."
679 int * index, *noe_index;
681 real ** d_r, **d, **dtot, **dtot2, **mean, **rms, **rmsc, *resnr;
682 real ** dtot1_3 = nullptr, **dtot1_6 = nullptr;
683 real rmsnow, meanmax, rmsmax, rmscmax;
685 t_noe_gr* noe_gr = nullptr;
686 t_noe** noe = nullptr;
688 gmx_bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
690 static int nlevels = 40;
691 static real scalemax = -1.0;
692 static gmx_bool bSumH = TRUE;
693 static gmx_bool bPBC = TRUE;
694 gmx_output_env_t* oenv;
697 { "-nlevels", FALSE, etINT, { &nlevels }, "Discretize RMS in this number of levels" },
698 { "-max", FALSE, etREAL, { &scalemax }, "Maximum level in matrices" },
699 { "-sumh", FALSE, etBOOL, { &bSumH }, "Average distance over equivalent hydrogens" },
704 "Use periodic boundary conditions when computing distances" }
707 { efTRX, "-f", nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
708 { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-equiv", "equiv", ffOPTRD },
709 { efXVG, nullptr, "distrmsd", ffWRITE }, { efXPM, "-rms", "rmsdist", ffOPTWR },
710 { efXPM, "-scl", "rmsscale", ffOPTWR }, { efXPM, "-mean", "rmsmean", ffOPTWR },
711 { efXPM, "-nmr3", "nmr3", ffOPTWR }, { efXPM, "-nmr6", "nmr6", ffOPTWR },
712 { efDAT, "-noe", "noe", ffOPTWR },
714 #define NFILE asize(fnm)
716 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
717 asize(desc), desc, 0, nullptr, &oenv))
722 bRMS = opt2bSet("-rms", NFILE, fnm);
723 bScale = opt2bSet("-scl", NFILE, fnm);
724 bMean = opt2bSet("-mean", NFILE, fnm);
725 bNOE = opt2bSet("-noe", NFILE, fnm);
726 bNMR3 = opt2bSet("-nmr3", NFILE, fnm);
727 bNMR6 = opt2bSet("-nmr6", NFILE, fnm);
728 bNMR = bNMR3 || bNMR6 || bNOE;
734 if (bNOE && scalemax < 0)
738 "WARNING: using -noe without -max "
739 "makes no sense, setting -max to %g\n\n",
743 /* get topology and index */
744 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box, FALSE);
750 atoms = &(top.atoms);
752 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
754 /* initialize arrays */
760 snew(dtot1_3, isize);
761 snew(dtot1_6, isize);
768 for (i = 0; (i < isize); i++)
771 snew(dtot[i], isize);
772 snew(dtot2[i], isize);
775 snew(dtot1_3[i], isize);
776 snew(dtot1_6[i], isize);
778 snew(mean[i], isize);
780 snew(rmsc[i], isize);
786 calc_dist(isize, index, x, ePBC, box, d_r);
789 /*open output files*/
790 fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS Deviation", "Time (ps)", "RMSD (nm)", oenv);
791 if (output_env_get_print_xvgr_codes(oenv))
793 fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
797 read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
802 calc_dist_tot(isize, index, x, ePBC, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
804 rmsnow = rms_diff(isize, d, d_r);
805 fprintf(fp, "%g %g\n", t, rmsnow);
807 } while (read_next_x(oenv, status, &t, x, box));
808 fprintf(stderr, "\n");
814 calc_rms(isize, teller, dtot, dtot2, rms, &rmsmax, rmsc, &rmscmax, mean, &meanmax);
815 fprintf(stderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
819 calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
833 /* make list of noe atom groups */
834 snew(noe_index, isize + 1);
836 gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE, fnm), atoms, isize, index, bSumH,
838 fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n", gnr, isize);
839 /* make half matrix of of noe-group distances from atom distances */
841 for (i = 0; i < gnr; i++)
845 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
857 write_xpm(opt2FILE("-rms", NFILE, fnm, "w"), 0, "RMS of distance", "RMS (nm)", "Atom Index",
858 "Atom Index", isize, isize, resnr, resnr, rms, 0.0, rmsmax, rlo, rhi, &nlevels);
863 write_xpm(opt2FILE("-scl", NFILE, fnm, "w"), 0, "Relative RMS", "RMS", "Atom Index",
864 "Atom Index", isize, isize, resnr, resnr, rmsc, 0.0, rmscmax, rlo, rhi, &nlevels);
869 write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean Distance", "Distance (nm)",
870 "Atom Index", "Atom Index", isize, isize, resnr, resnr, mean, 0.0, meanmax, rlo,
876 write_xpm(opt2FILE("-nmr3", NFILE, fnm, "w"), 0, "1/r^3 averaged distances",
877 "Distance (nm)", "Atom Index", "Atom Index", isize, isize, resnr, resnr, dtot1_3,
878 0.0, max1_3, rlo, rhi, &nlevels);
882 write_xpm(opt2FILE("-nmr6", NFILE, fnm, "w"), 0, "1/r^6 averaged distances",
883 "Distance (nm)", "Atom Index", "Atom Index", isize, isize, resnr, resnr, dtot1_6,
884 0.0, max1_6, rlo, rhi, &nlevels);
889 write_noe(opt2FILE("-noe", NFILE, fnm, "w"), gnr, noe, noe_gr, scalemax);
892 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);