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, 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.
41 #include "gromacs/commandline/pargs.h"
42 #include "gromacs/fileio/matio.h"
43 #include "gromacs/fileio/strdb.h"
44 #include "gromacs/fileio/tpxio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/gmx_ana.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/legacyheaders/names.h"
50 #include "gromacs/legacyheaders/typedefs.h"
51 #include "gromacs/legacyheaders/viewit.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/pbcutil/pbc.h"
54 #include "gromacs/topology/index.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
59 static void calc_dist(int nind, atom_id index[], rvec x[], int ePBC, matrix box,
68 set_pbc(&pbc, ePBC, box);
69 for (i = 0; (i < nind-1); i++)
72 for (j = i+1; (j < nind); j++)
74 pbc_dx(&pbc, xi, x[index[j]], dx);
76 d[i][j] = sqrt(temp2);
81 static void calc_dist_tot(int nind, atom_id index[], rvec x[],
83 real **d, real **dtot, real **dtot2,
84 gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
88 real temp, temp2, temp1_3;
92 set_pbc(&pbc, ePBC, box);
93 for (i = 0; (i < nind-1); i++)
96 for (j = i+1; (j < nind); j++)
98 pbc_dx(&pbc, xi, x[index[j]], dx);
103 dtot2[i][j] += temp2;
106 temp1_3 = 1.0/(temp*temp2);
107 dtot1_3[i][j] += temp1_3;
108 dtot1_6[i][j] += temp1_3*temp1_3;
114 static void calc_nmr(int nind, int nframes, real **dtot1_3, real **dtot1_6,
115 real *max1_3, real *max1_6)
118 real temp1_3, temp1_6;
120 for (i = 0; (i < nind-1); i++)
122 for (j = i+1; (j < nind); j++)
124 temp1_3 = pow(dtot1_3[i][j]/nframes, -1.0/3.0);
125 temp1_6 = pow(dtot1_6[i][j]/nframes, -1.0/6.0);
126 if (temp1_3 > *max1_3)
130 if (temp1_6 > *max1_6)
134 dtot1_3[i][j] = temp1_3;
135 dtot1_6[i][j] = temp1_6;
136 dtot1_3[j][i] = temp1_3;
137 dtot1_6[j][i] = temp1_6;
142 static char Hnum[] = "123";
167 static int read_equiv(const char *eq_fn, t_equiv ***equivptr)
170 char line[STRLEN], resname[10], atomname[10], *lp;
171 int neq, na, n, resnr;
174 fp = gmx_ffopen(eq_fn, "r");
177 while (get_a_line(fp, line, STRLEN))
180 /* this is not efficient, but I'm lazy */
181 srenew(equiv, neq+1);
184 if (sscanf(lp, "%s %n", atomname, &n) == 1)
188 equiv[neq][0].nname = gmx_strdup(atomname);
189 while (sscanf(lp, "%d %s %s %n", &resnr, resname, atomname, &n) == 3)
191 /* this is not efficient, but I'm lazy (again) */
192 srenew(equiv[neq], na+1);
193 equiv[neq][na].rnr = resnr-1;
194 equiv[neq][na].rname = gmx_strdup(resname);
195 equiv[neq][na].aname = gmx_strdup(atomname);
198 equiv[neq][na].nname = NULL;
204 /* make empty element as flag for end of array */
205 srenew(equiv[neq], na+1);
206 equiv[neq][na].rnr = NOTSET;
207 equiv[neq][na].rname = NULL;
208 equiv[neq][na].aname = NULL;
220 static void dump_equiv(FILE *out, int neq, t_equiv **equiv)
224 fprintf(out, "Dumping equivalent list\n");
225 for (i = 0; i < neq; i++)
227 fprintf(out, "%s", equiv[i][0].nname);
228 for (j = 0; equiv[i][j].rnr != NOTSET; j++)
230 fprintf(out, " %d %s %s",
231 equiv[i][j].rnr, equiv[i][j].rname, equiv[i][j].aname);
237 static gmx_bool is_equiv(int neq, t_equiv **equiv, char **nname,
238 int rnr1, char *rname1, char *aname1,
239 int rnr2, char *rname2, char *aname2)
245 /* we can terminate each loop when bFound is true! */
246 for (i = 0; i < neq && !bFound; i++)
248 /* find first atom */
249 for (j = 0; equiv[i][j].rnr != NOTSET && !bFound; j++)
251 bFound = ( equiv[i][j].rnr == rnr1 &&
252 strcmp(equiv[i][j].rname, rname1) == 0 &&
253 strcmp(equiv[i][j].aname, aname1) == 0 );
257 /* find second atom */
259 for (j = 0; equiv[i][j].rnr != NOTSET && !bFound; j++)
261 bFound = ( equiv[i][j].rnr == rnr2 &&
262 strcmp(equiv[i][j].rname, rname2) == 0 &&
263 strcmp(equiv[i][j].aname, aname2) == 0 );
269 *nname = gmx_strdup(equiv[i-1][0].nname);
275 static int analyze_noe_equivalent(const char *eq_fn,
276 t_atoms *atoms, int isize, atom_id *index,
278 atom_id *noe_index, t_noe_gr *noe_gr)
281 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
282 char *anmi, *anmj, **nnm;
283 gmx_bool bMatch, bEquiv;
291 neq = read_equiv(eq_fn, &equiv);
294 dump_equiv(debug, neq, equiv);
304 for (i = 0; i < isize; i++)
306 if (equiv && i < isize-1)
308 /* check explicit list of equivalent atoms */
312 rnri = atoms->atom[index[i]].resind;
313 rnrj = atoms->atom[index[j]].resind;
315 is_equiv(neq, equiv, &nnm[i],
316 rnri, *atoms->resinfo[rnri].name, *atoms->atomname[index[i]],
317 rnrj, *atoms->resinfo[rnrj].name, *atoms->atomname[index[j]]);
318 if (nnm[i] && bEquiv)
320 nnm[j] = gmx_strdup(nnm[i]);
324 /* set index for matching atom */
325 noe_index[i] = groupnr;
326 /* skip matching atom */
330 while (bEquiv && i < isize-1);
338 /* look for triplets of consecutive atoms with name XX?,
339 X are any number of letters or digits and ? goes from 1 to 3
340 This is supposed to cover all CH3 groups and the like */
341 anmi = *atoms->atomname[index[i]];
342 anmil = strlen(anmi);
343 bMatch = i <= isize-3 && anmi[anmil-1] == '1';
346 for (j = 1; j < 3; j++)
348 anmj = *atoms->atomname[index[i+j]];
349 anmjl = strlen(anmj);
350 bMatch = bMatch && ( anmil == anmjl && anmj[anmjl-1] == Hnum[j] &&
351 strncmp(anmi, anmj, anmil-1) == 0 );
354 /* set index for this atom */
355 noe_index[i] = groupnr;
358 /* set index for next two matching atoms */
359 for (j = 1; j < 3; j++)
361 noe_index[i+j] = groupnr;
363 /* skip matching atoms */
372 /* make index without looking for equivalent atoms */
373 for (i = 0; i < isize; i++)
379 noe_index[isize] = groupnr;
384 for (i = 0; i < isize; i++)
386 rnri = atoms->atom[index[i]].resind;
387 fprintf(debug, "%s %s %d -> %s\n", *atoms->atomname[index[i]],
388 *atoms->resinfo[rnri].name, rnri, nnm[i] ? nnm[i] : "");
392 for (i = 0; i < isize; i++)
395 if (!noe_gr[gi].aname)
398 noe_gr[gi].anr = index[i];
401 noe_gr[gi].aname = gmx_strdup(nnm[i]);
405 noe_gr[gi].aname = gmx_strdup(*atoms->atomname[index[i]]);
406 if (noe_index[i] == noe_index[i+1])
408 noe_gr[gi].aname[strlen(noe_gr[gi].aname)-1] = '*';
411 noe_gr[gi].rnr = atoms->atom[index[i]].resind;
412 noe_gr[gi].rname = gmx_strdup(*atoms->resinfo[noe_gr[gi].rnr].name);
413 /* dump group definitions */
416 fprintf(debug, "%d %d %d %d %s %s %d\n", i, gi,
417 noe_gr[gi].ianr, noe_gr[gi].anr, noe_gr[gi].aname,
418 noe_gr[gi].rname, noe_gr[gi].rnr);
422 for (i = 0; i < isize; i++)
431 /* #define NSCALE 3 */
432 /* static char *noe_scale[NSCALE+1] = { */
433 /* "strong", "medium", "weak", "none" */
437 static char *noe2scale(real r3, real r6, real rmax)
440 static char buf[NSCALE+1];
442 /* r goes from 0 to rmax
443 NSCALE*r/rmax goes from 0 to NSCALE
444 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
445 s3 = NSCALE - min(NSCALE, (int)(NSCALE*r3/rmax));
446 s6 = NSCALE - min(NSCALE, (int)(NSCALE*r6/rmax));
448 for (i = 0; i < s3; i++)
461 static void calc_noe(int isize, atom_id *noe_index,
462 real **dtot1_3, real **dtot1_6, int gnr, t_noe **noe)
466 /* make half matrix of noe-group distances from atom distances */
467 for (i = 0; i < isize; i++)
470 for (j = i; j < isize; j++)
474 noe[gi][gj].i_3 += pow(dtot1_3[i][j], -3);
475 noe[gi][gj].i_6 += pow(dtot1_6[i][j], -6);
480 for (i = 0; i < gnr; i++)
482 for (j = i+1; j < gnr; j++)
484 noe[i][j].r_3 = pow(noe[i][j].i_3/noe[i][j].nr, -1.0/3.0);
485 noe[i][j].r_6 = pow(noe[i][j].i_6/noe[i][j].nr, -1.0/6.0);
486 noe[j][i] = noe[i][j];
491 static void write_noe(FILE *fp, int gnr, t_noe **noe, t_noe_gr *noe_gr, real rmax)
494 real r3, r6, min3, min6;
495 char buf[10], b3[10], b6[10];
500 ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
501 "ianr", "anr", "anm", "rnm", "rnr", "ianr", "anr", "anm", "rnm", "rnr",
502 "1/r^3", "1/r^6", "intnsty", "Dr", "Da", "scale");
503 for (i = 0; i < gnr; i++)
506 for (j = i+1; j < gnr; j++)
511 min3 = min(r3, min3);
512 min6 = min(r6, min6);
513 if (r3 < rmax || r6 < rmax)
515 if (grj.rnr == gri.rnr)
517 sprintf(buf, "%2d", grj.anr-gri.anr);
525 sprintf(b3, "%-5.3f", r3);
533 sprintf(b6, "%-5.3f", r6);
540 "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
541 gri.ianr+1, gri.anr+1, gri.aname, gri.rname, gri.rnr+1,
542 grj.ianr+1, grj.anr+1, grj.aname, grj.rname, grj.rnr+1,
543 b3, b6, (int)(noe[i][j].i_6+0.5), grj.rnr-gri.rnr, buf,
544 noe2scale(r3, r6, rmax));
548 #define MINI ((i == 3) ? min3 : min6)
549 for (i = 3; i <= 6; i += 3)
553 fprintf(stdout, "NOTE: no 1/r^%d averaged distances found below %g, "
554 "smallest was %g\n", i, rmax, MINI );
558 fprintf(stdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI );
564 static void calc_rms(int nind, int nframes,
565 real **dtot, real **dtot2,
566 real **rmsmat, real *rmsmax,
567 real **rmscmat, real *rmscmax,
568 real **meanmat, real *meanmax)
571 real mean, mean2, rms, rmsc;
572 /* N.B. dtot and dtot2 contain the total distance and the total squared
573 * distance respectively, BUT they return RMS and the scaled RMS resp.
580 for (i = 0; (i < nind-1); i++)
582 for (j = i+1; (j < nind); j++)
584 mean = dtot[i][j]/nframes;
585 mean2 = dtot2[i][j]/nframes;
586 rms = sqrt(max(0, mean2-mean*mean));
600 meanmat[i][j] = meanmat[j][i] = mean;
601 rmsmat[i][j] = rmsmat[j][i] = rms;
602 rmscmat[i][j] = rmscmat[j][i] = rmsc;
607 real rms_diff(int natom, real **d, real **d_r)
613 for (i = 0; (i < natom-1); i++)
615 for (j = i+1; (j < natom); j++)
617 r = d[i][j]-d_r[i][j];
621 r2 /= (natom*(natom-1))/2;
626 int gmx_rmsdist(int argc, char *argv[])
628 const char *desc[] = {
629 "[THISMODULE] computes the root mean square deviation of atom distances,",
630 "which has the advantage that no fit is needed like in standard RMS",
631 "deviation as computed by [gmx-rms].",
632 "The reference structure is taken from the structure file.",
633 "The RMSD at time t is calculated as the RMS",
634 "of the differences in distance between atom-pairs in the reference",
635 "structure and the structure at time t.[PAR]",
636 "[THISMODULE] can also produce matrices of the rms distances, rms distances",
637 "scaled with the mean distance and the mean distances and matrices with",
638 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
639 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
640 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
641 "can be generated, by default averaging over equivalent hydrogens",
642 "(all triplets of hydrogens named *[123]). Additionally a list of",
643 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
644 "a set of equivalent atoms specified as residue number and name and",
645 "atom name; e.g.:[PAR]",
646 "[TT]HB* 3 SER HB1 3 SER HB2[tt][PAR]",
647 "Residue and atom names must exactly match those in the structure",
648 "file, including case. Specifying non-sequential atoms is undefined."
652 int natom, i, j, teller, gi, gj;
664 atom_id *index, *noe_index;
666 real **d_r, **d, **dtot, **dtot2, **mean, **rms, **rmsc, *resnr;
667 real **dtot1_3 = NULL, **dtot1_6 = NULL;
668 real rmsnow, meanmax, rmsmax, rmscmax;
670 t_noe_gr *noe_gr = NULL;
674 gmx_bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
676 static int nlevels = 40;
677 static real scalemax = -1.0;
678 static gmx_bool bSumH = TRUE;
679 static gmx_bool bPBC = TRUE;
683 { "-nlevels", FALSE, etINT, {&nlevels},
684 "Discretize RMS in this number of levels" },
685 { "-max", FALSE, etREAL, {&scalemax},
686 "Maximum level in matrices" },
687 { "-sumh", FALSE, etBOOL, {&bSumH},
688 "Average distance over equivalent hydrogens" },
689 { "-pbc", FALSE, etBOOL, {&bPBC},
690 "Use periodic boundary conditions when computing distances" }
693 { efTRX, "-f", NULL, ffREAD },
694 { efTPS, NULL, NULL, ffREAD },
695 { efNDX, NULL, NULL, ffOPTRD },
696 { efDAT, "-equiv", "equiv", ffOPTRD },
697 { efXVG, NULL, "distrmsd", ffWRITE },
698 { efXPM, "-rms", "rmsdist", ffOPTWR },
699 { efXPM, "-scl", "rmsscale", ffOPTWR },
700 { efXPM, "-mean", "rmsmean", ffOPTWR },
701 { efXPM, "-nmr3", "nmr3", ffOPTWR },
702 { efXPM, "-nmr6", "nmr6", ffOPTWR },
703 { efDAT, "-noe", "noe", ffOPTWR },
705 #define NFILE asize(fnm)
707 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
708 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
713 bRMS = opt2bSet("-rms", NFILE, fnm);
714 bScale = opt2bSet("-scl", NFILE, fnm);
715 bMean = opt2bSet("-mean", NFILE, fnm);
716 bNOE = opt2bSet("-noe", NFILE, fnm);
717 bNMR3 = opt2bSet("-nmr3", NFILE, fnm);
718 bNMR6 = opt2bSet("-nmr6", NFILE, fnm);
719 bNMR = bNMR3 || bNMR6 || bNOE;
725 if (bNOE && scalemax < 0)
728 fprintf(stderr, "WARNING: using -noe without -max "
729 "makes no sense, setting -max to %g\n\n", scalemax);
732 /* get topology and index */
733 read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &x, NULL, box, FALSE);
739 atoms = &(top.atoms);
741 get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
743 /* initialize arrays */
749 snew(dtot1_3, isize);
750 snew(dtot1_6, isize);
757 for (i = 0; (i < isize); i++)
760 snew(dtot[i], isize);
761 snew(dtot2[i], isize);
764 snew(dtot1_3[i], isize);
765 snew(dtot1_6[i], isize);
767 snew(mean[i], isize);
769 snew(rmsc[i], isize);
775 calc_dist(isize, index, x, ePBC, box, d_r);
778 /*open output files*/
779 fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS Deviation", "Time (ps)", "RMSD (nm)",
781 if (output_env_get_print_xvgr_codes(oenv))
783 fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
787 natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
792 calc_dist_tot(isize, index, x, ePBC, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
794 rmsnow = rms_diff(isize, d, d_r);
795 fprintf(fp, "%g %g\n", t, rmsnow);
798 while (read_next_x(oenv, status, &t, x, box));
799 fprintf(stderr, "\n");
805 calc_rms(isize, teller, dtot, dtot2, mean, &meanmax, rms, &rmsmax, rmsc, &rmscmax);
806 fprintf(stderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
810 calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
824 /* make list of noe atom groups */
825 snew(noe_index, isize+1);
827 gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE, fnm),
828 atoms, isize, index, bSumH, noe_index, noe_gr);
829 fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n",
831 /* make half matrix of of noe-group distances from atom distances */
833 for (i = 0; i < gnr; i++)
837 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
840 rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
841 rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;
845 write_xpm(opt2FILE("-rms", NFILE, fnm, "w"), 0,
846 "RMS of distance", "RMS (nm)", "Atom Index", "Atom Index",
847 isize, isize, resnr, resnr, rms, 0.0, rmsmax, rlo, rhi, &nlevels);
852 write_xpm(opt2FILE("-scl", NFILE, fnm, "w"), 0,
853 "Relative RMS", "RMS", "Atom Index", "Atom Index",
854 isize, isize, resnr, resnr, rmsc, 0.0, rmscmax, rlo, rhi, &nlevels);
859 write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0,
860 "Mean Distance", "Distance (nm)", "Atom Index", "Atom Index",
861 isize, isize, resnr, resnr, mean, 0.0, meanmax, rlo, rhi, &nlevels);
866 write_xpm(opt2FILE("-nmr3", NFILE, fnm, "w"), 0, "1/r^3 averaged distances",
867 "Distance (nm)", "Atom Index", "Atom Index",
868 isize, isize, resnr, resnr, dtot1_3, 0.0, max1_3, rlo, rhi, &nlevels);
872 write_xpm(opt2FILE("-nmr6", NFILE, fnm, "w"), 0, "1/r^6 averaged distances",
873 "Distance (nm)", "Atom Index", "Atom Index",
874 isize, isize, resnr, resnr, dtot1_6, 0.0, max1_6, rlo, rhi, &nlevels);
879 write_noe(opt2FILE("-noe", NFILE, fnm, "w"), gnr, noe, noe_gr, scalemax);
882 do_view(oenv, ftp2fn(efXVG, NFILE, fnm), NULL);