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 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.
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.
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.
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.
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.
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.
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/gmx_ana.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/topology/index.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/utility/arraysize.h"
59 #include "gromacs/utility/cstringutil.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/futil.h"
62 #include "gromacs/utility/smalloc.h"
63 #include "gromacs/utility/strdb.h"
66 static void calc_dist(int nind, const int index[], const rvec x[], PbcType pbcType, matrix box, real** d)
73 set_pbc(&pbc, pbcType, box);
74 for (i = 0; (i < nind - 1); i++)
76 const real* xi = x[index[i]];
77 for (j = i + 1; (j < nind); j++)
79 pbc_dx(&pbc, xi, x[index[j]], dx);
81 d[i][j] = std::sqrt(temp2);
86 static void calc_dist_tot(int nind,
100 real temp, temp2, temp1_3;
104 set_pbc(&pbc, pbcType, box);
105 for (i = 0; (i < nind - 1); i++)
108 for (j = i + 1; (j < nind); j++)
110 pbc_dx(&pbc, xi, x[index[j]], dx);
112 temp = std::sqrt(temp2);
115 dtot2[i][j] += temp2;
118 temp1_3 = 1.0 / (temp * temp2);
119 dtot1_3[i][j] += temp1_3;
120 dtot1_6[i][j] += temp1_3 * temp1_3;
126 static void calc_nmr(int nind, int nframes, real** dtot1_3, real** dtot1_6, real* max1_3, real* max1_6)
129 real temp1_3, temp1_6;
131 for (i = 0; (i < nind - 1); i++)
133 for (j = i + 1; (j < nind); j++)
135 temp1_3 = gmx::invcbrt(dtot1_3[i][j] / static_cast<real>(nframes));
136 temp1_6 = gmx::invsixthroot(dtot1_6[i][j] / static_cast<real>(nframes));
137 if (temp1_3 > *max1_3)
141 if (temp1_6 > *max1_6)
145 dtot1_3[i][j] = temp1_3;
146 dtot1_6[i][j] = temp1_6;
147 dtot1_3[j][i] = temp1_3;
148 dtot1_6[j][i] = temp1_6;
153 static char Hnum[] = "123";
182 static int read_equiv(const char* eq_fn, t_equiv*** equivptr)
185 char line[STRLEN], resname[10], atomname[10], *lp;
186 int neq, na, n, resnr;
189 fp = gmx_ffopen(eq_fn, "r");
192 while (get_a_line(fp, line, STRLEN))
195 /* this is not efficient, but I'm lazy */
196 srenew(equiv, neq + 1);
197 equiv[neq] = nullptr;
199 if (sscanf(lp, "%s %n", atomname, &n) == 1)
203 equiv[neq][0].nname = gmx_strdup(atomname);
204 while (sscanf(lp, "%d %s %s %n", &resnr, resname, atomname, &n) == 3)
206 /* this is not efficient, but I'm lazy (again) */
207 srenew(equiv[neq], na + 1);
208 equiv[neq][na].set = true;
209 equiv[neq][na].rnr = resnr - 1;
210 equiv[neq][na].rname = gmx_strdup(resname);
211 equiv[neq][na].aname = gmx_strdup(atomname);
214 equiv[neq][na].nname = nullptr;
220 /* make empty element as flag for end of array */
221 srenew(equiv[neq], na + 1);
222 equiv[neq][na].set = false;
223 equiv[neq][na].rnr = 0;
224 equiv[neq][na].rname = nullptr;
225 equiv[neq][na].aname = nullptr;
237 static void dump_equiv(FILE* out, int neq, t_equiv** equiv)
241 fprintf(out, "Dumping equivalent list\n");
242 for (i = 0; i < neq; i++)
244 fprintf(out, "%s", equiv[i][0].nname);
245 for (j = 0; equiv[i][j].set; j++)
247 fprintf(out, " %d %s %s", equiv[i][j].rnr, equiv[i][j].rname, equiv[i][j].aname);
254 is_equiv(int neq, t_equiv** equiv, char** nname, int rnr1, char* rname1, char* aname1, int rnr2, char* rname2, char* aname2)
260 /* we can terminate each loop when bFound is true! */
261 for (i = 0; i < neq && !bFound; i++)
263 /* find first atom */
264 for (j = 0; equiv[i][j].set && !bFound; j++)
266 bFound = (equiv[i][j].rnr == rnr1 && std::strcmp(equiv[i][j].rname, rname1) == 0
267 && std::strcmp(equiv[i][j].aname, aname1) == 0);
271 /* find second atom */
273 for (j = 0; equiv[i][j].set && !bFound; j++)
275 bFound = (equiv[i][j].rnr == rnr2 && std::strcmp(equiv[i][j].rname, rname2) == 0
276 && std::strcmp(equiv[i][j].aname, aname2) == 0);
282 *nname = gmx_strdup(equiv[i - 1][0].nname);
288 static int analyze_noe_equivalent(const char* eq_fn,
289 const t_atoms* atoms,
296 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
297 char * anmi, *anmj, **nnm;
298 gmx_bool bMatch, bEquiv;
306 neq = read_equiv(eq_fn, &equiv);
309 dump_equiv(debug, neq, equiv);
319 for (i = 0; i < isize; i++)
321 if (equiv && i < isize - 1)
323 /* check explicit list of equivalent atoms */
327 rnri = atoms->atom[index[i]].resind;
328 rnrj = atoms->atom[index[j]].resind;
329 bEquiv = is_equiv(neq, equiv, &nnm[i], rnri, *atoms->resinfo[rnri].name,
330 *atoms->atomname[index[i]], rnrj, *atoms->resinfo[rnrj].name,
331 *atoms->atomname[index[j]]);
332 if (nnm[i] && bEquiv)
334 nnm[j] = gmx_strdup(nnm[i]);
338 /* set index for matching atom */
339 noe_index[i] = groupnr;
340 /* skip matching atom */
343 } while (bEquiv && i < isize - 1);
351 /* look for triplets of consecutive atoms with name XX?,
352 X are any number of letters or digits and ? goes from 1 to 3
353 This is supposed to cover all CH3 groups and the like */
354 anmi = *atoms->atomname[index[i]];
355 anmil = std::strlen(anmi);
356 bMatch = i <= isize - 3 && anmi[anmil - 1] == '1';
359 for (j = 1; j < 3; j++)
361 anmj = *atoms->atomname[index[i + j]];
362 anmjl = std::strlen(anmj);
364 && (anmil == anmjl && anmj[anmjl - 1] == Hnum[j]
365 && std::strncmp(anmi, anmj, anmil - 1) == 0);
368 /* set index for this atom */
369 noe_index[i] = groupnr;
372 /* set index for next two matching atoms */
373 for (j = 1; j < 3; j++)
375 noe_index[i + j] = groupnr;
377 /* skip matching atoms */
386 /* make index without looking for equivalent atoms */
387 for (i = 0; i < isize; i++)
393 noe_index[isize] = groupnr;
398 for (i = 0; i < isize; i++)
400 rnri = atoms->atom[index[i]].resind;
401 fprintf(debug, "%s %s %d -> %s\n", *atoms->atomname[index[i]],
402 *atoms->resinfo[rnri].name, rnri, nnm[i] ? nnm[i] : "");
406 for (i = 0; i < isize; i++)
409 if (!noe_gr[gi].aname)
412 noe_gr[gi].anr = index[i];
415 noe_gr[gi].aname = gmx_strdup(nnm[i]);
419 noe_gr[gi].aname = gmx_strdup(*atoms->atomname[index[i]]);
420 if (noe_index[i] == noe_index[i + 1])
422 noe_gr[gi].aname[std::strlen(noe_gr[gi].aname) - 1] = '*';
425 noe_gr[gi].rnr = atoms->atom[index[i]].resind;
426 noe_gr[gi].rname = gmx_strdup(*atoms->resinfo[noe_gr[gi].rnr].name);
427 /* dump group definitions */
430 fprintf(debug, "%d %d %d %d %s %s %d\n", i, gi, noe_gr[gi].ianr, noe_gr[gi].anr,
431 noe_gr[gi].aname, noe_gr[gi].rname, noe_gr[gi].rnr);
435 for (i = 0; i < isize; i++)
444 /* #define NSCALE 3 */
445 /* static char *noe_scale[NSCALE+1] = { */
446 /* "strong", "medium", "weak", "none" */
450 static char* noe2scale(real r3, real r6, real rmax)
453 static char buf[NSCALE + 1];
455 /* r goes from 0 to rmax
456 NSCALE*r/rmax goes from 0 to NSCALE
457 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
458 s3 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE * r3 / rmax));
459 s6 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE * r6 / rmax));
461 for (i = 0; i < s3; i++)
474 static void calc_noe(int isize, const int* noe_index, real** dtot1_3, real** dtot1_6, int gnr, t_noe** noe)
478 /* make half matrix of noe-group distances from atom distances */
479 for (i = 0; i < isize; i++)
482 for (j = i; j < isize; j++)
486 noe[gi][gj].i_3 += 1.0 / gmx::power3(dtot1_3[i][j]);
487 noe[gi][gj].i_6 += 1.0 / gmx::power6(dtot1_6[i][j]);
492 for (i = 0; i < gnr; i++)
494 for (j = i + 1; j < gnr; j++)
496 noe[i][j].r_3 = gmx::invcbrt(noe[i][j].i_3 / static_cast<real>(noe[i][j].nr));
497 noe[i][j].r_6 = gmx::invsixthroot(noe[i][j].i_6 / static_cast<real>(noe[i][j].nr));
498 noe[j][i] = noe[i][j];
503 static void write_noe(FILE* fp, int gnr, t_noe** noe, t_noe_gr* noe_gr, real rmax)
506 real r3, r6, min3, min6;
507 char buf[10], b3[10], b6[10];
511 fprintf(fp, ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n", "ianr", "anr",
512 "anm", "rnm", "rnr", "ianr", "anr", "anm", "rnm", "rnr", "1/r^3", "1/r^6", "intnsty",
513 "Dr", "Da", "scale");
514 for (i = 0; i < gnr; i++)
517 for (j = i + 1; j < gnr; j++)
522 min3 = std::min(r3, min3);
523 min6 = std::min(r6, min6);
524 if (r3 < rmax || r6 < rmax)
526 if (grj.rnr == gri.rnr)
528 sprintf(buf, "%2d", grj.anr - gri.anr);
536 sprintf(b3, "%-5.3f", r3);
540 std::strcpy(b3, "-");
544 sprintf(b6, "%-5.3f", r6);
548 std::strcpy(b6, "-");
550 fprintf(fp, "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
551 gri.ianr + 1, gri.anr + 1, gri.aname, gri.rname, gri.rnr + 1, grj.ianr + 1,
552 grj.anr + 1, grj.aname, grj.rname, grj.rnr + 1, b3, b6,
553 gmx::roundToInt(noe[i][j].i_6), grj.rnr - gri.rnr, buf, noe2scale(r3, r6, rmax));
557 #define MINI ((i == 3) ? min3 : min6)
558 for (i = 3; i <= 6; i += 3)
563 "NOTE: no 1/r^%d averaged distances found below %g, "
569 fprintf(stdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI);
575 static void calc_rms(int nind,
587 real mean, mean2, rms, rmsc;
588 /* N.B. dtot and dtot2 contain the total distance and the total squared
589 * distance respectively, BUT they return RMS and the scaled RMS resp.
596 for (i = 0; (i < nind - 1); i++)
598 for (j = i + 1; (j < nind); j++)
600 mean = dtot[i][j] / static_cast<real>(nframes);
601 mean2 = dtot2[i][j] / static_cast<real>(nframes);
602 rms = std::sqrt(std::max(0.0_real, mean2 - mean * mean));
616 meanmat[i][j] = meanmat[j][i] = mean;
617 rmsmat[i][j] = rmsmat[j][i] = rms;
618 rmscmat[i][j] = rmscmat[j][i] = rmsc;
623 static real rms_diff(int natom, real** d, real** d_r)
629 for (i = 0; (i < natom - 1); i++)
631 for (j = i + 1; (j < natom); j++)
633 r = d[i][j] - d_r[i][j];
637 r2 /= gmx::exactDiv(natom * (natom - 1), 2);
639 return std::sqrt(r2);
642 int gmx_rmsdist(int argc, char* argv[])
644 const char* desc[] = {
645 "[THISMODULE] computes the root mean square deviation of atom distances,",
646 "which has the advantage that no fit is needed like in standard RMS",
647 "deviation as computed by [gmx-rms].",
648 "The reference structure is taken from the structure file.",
649 "The RMSD at time t is calculated as the RMS",
650 "of the differences in distance between atom-pairs in the reference",
651 "structure and the structure at time t.[PAR]",
652 "[THISMODULE] can also produce matrices of the rms distances, rms distances",
653 "scaled with the mean distance and the mean distances and matrices with",
654 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
655 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
656 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
657 "can be generated, by default averaging over equivalent hydrogens",
658 "(all triplets of hydrogens named \\*[123]). Additionally a list of",
659 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
660 "a set of equivalent atoms specified as residue number and name and",
661 "atom name; e.g.:[PAR]",
662 "[TT]HB* 3 SER HB1 3 SER HB2[tt][PAR]",
663 "Residue and atom names must exactly match those in the structure",
664 "file, including case. Specifying non-sequential atoms is undefined."
680 int * index, *noe_index;
682 real ** d_r, **d, **dtot, **dtot2, **mean, **rms, **rmsc, *resnr;
683 real ** dtot1_3 = nullptr, **dtot1_6 = nullptr;
684 real rmsnow, meanmax, rmsmax, rmscmax;
686 t_noe_gr* noe_gr = nullptr;
687 t_noe** noe = nullptr;
689 gmx_bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
691 static int nlevels = 40;
692 static real scalemax = -1.0;
693 static gmx_bool bSumH = TRUE;
694 static gmx_bool bPBC = TRUE;
695 gmx_output_env_t* oenv;
698 { "-nlevels", FALSE, etINT, { &nlevels }, "Discretize RMS in this number of levels" },
699 { "-max", FALSE, etREAL, { &scalemax }, "Maximum level in matrices" },
700 { "-sumh", FALSE, etBOOL, { &bSumH }, "Average distance over equivalent hydrogens" },
705 "Use periodic boundary conditions when computing distances" }
708 { efTRX, "-f", nullptr, ffREAD }, { efTPS, nullptr, nullptr, ffREAD },
709 { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-equiv", "equiv", ffOPTRD },
710 { efXVG, nullptr, "distrmsd", ffWRITE }, { efXPM, "-rms", "rmsdist", ffOPTWR },
711 { efXPM, "-scl", "rmsscale", ffOPTWR }, { efXPM, "-mean", "rmsmean", ffOPTWR },
712 { efXPM, "-nmr3", "nmr3", ffOPTWR }, { efXPM, "-nmr6", "nmr6", ffOPTWR },
713 { efDAT, "-noe", "noe", ffOPTWR },
715 #define NFILE asize(fnm)
717 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
718 asize(desc), desc, 0, nullptr, &oenv))
723 bRMS = opt2bSet("-rms", NFILE, fnm);
724 bScale = opt2bSet("-scl", NFILE, fnm);
725 bMean = opt2bSet("-mean", NFILE, fnm);
726 bNOE = opt2bSet("-noe", NFILE, fnm);
727 bNMR3 = opt2bSet("-nmr3", NFILE, fnm);
728 bNMR6 = opt2bSet("-nmr6", NFILE, fnm);
729 bNMR = bNMR3 || bNMR6 || bNOE;
735 if (bNOE && scalemax < 0)
739 "WARNING: using -noe without -max "
740 "makes no sense, setting -max to %g\n\n",
744 /* get topology and index */
745 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &x, nullptr, box, FALSE);
749 pbcType = PbcType::No;
751 atoms = &(top.atoms);
753 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
755 /* initialize arrays */
761 snew(dtot1_3, isize);
762 snew(dtot1_6, isize);
769 for (i = 0; (i < isize); i++)
772 snew(dtot[i], isize);
773 snew(dtot2[i], isize);
776 snew(dtot1_3[i], isize);
777 snew(dtot1_6[i], isize);
779 snew(mean[i], isize);
781 snew(rmsc[i], isize);
787 calc_dist(isize, index, x, pbcType, box, d_r);
790 /*open output files*/
791 fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS Deviation", "Time (ps)", "RMSD (nm)", oenv);
792 if (output_env_get_print_xvgr_codes(oenv))
794 fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
798 read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
803 calc_dist_tot(isize, index, x, pbcType, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
805 rmsnow = rms_diff(isize, d, d_r);
806 fprintf(fp, "%g %g\n", t, rmsnow);
808 } while (read_next_x(oenv, status, &t, x, box));
809 fprintf(stderr, "\n");
815 calc_rms(isize, teller, dtot, dtot2, rms, &rmsmax, rmsc, &rmscmax, mean, &meanmax);
816 fprintf(stderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
820 calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
834 /* make list of noe atom groups */
835 snew(noe_index, isize + 1);
837 gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE, fnm), atoms, isize, index, bSumH,
839 fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n", gnr, isize);
840 /* make half matrix of of noe-group distances from atom distances */
842 for (i = 0; i < gnr; i++)
846 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
858 write_xpm(opt2FILE("-rms", NFILE, fnm, "w"), 0, "RMS of distance", "RMS (nm)", "Atom Index",
859 "Atom Index", isize, isize, resnr, resnr, rms, 0.0, rmsmax, rlo, rhi, &nlevels);
864 write_xpm(opt2FILE("-scl", NFILE, fnm, "w"), 0, "Relative RMS", "RMS", "Atom Index",
865 "Atom Index", isize, isize, resnr, resnr, rmsc, 0.0, rmscmax, rlo, rhi, &nlevels);
870 write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean Distance", "Distance (nm)",
871 "Atom Index", "Atom Index", isize, isize, resnr, resnr, mean, 0.0, meanmax, rlo,
877 write_xpm(opt2FILE("-nmr3", NFILE, fnm, "w"), 0, "1/r^3 averaged distances",
878 "Distance (nm)", "Atom Index", "Atom Index", isize, isize, resnr, resnr, dtot1_3,
879 0.0, max1_3, rlo, rhi, &nlevels);
883 write_xpm(opt2FILE("-nmr6", NFILE, fnm, "w"), 0, "1/r^6 averaged distances",
884 "Distance (nm)", "Atom Index", "Atom Index", isize, isize, resnr, resnr, dtot1_6,
885 0.0, max1_6, rlo, rhi, &nlevels);
890 write_noe(opt2FILE("-noe", NFILE, fnm, "w"), gnr, noe, noe_gr, scalemax);
893 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);