src/gromacs/gmxana/gmx_rmsdist.cpp

   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, 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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstring>
  41
  42 #include <algorithm>
  43
  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"
  63
  64
  65 static void calc_dist(int nind, int index[], const rvec x[], int ePBC, matrix box,
  66                       real **d)
  67 {
  68     int      i, j;
  69     rvec     dx;
  70     real     temp2;
  71     t_pbc    pbc;
  72
  73     set_pbc(&pbc, ePBC, box);
  74     for (i = 0; (i < nind-1); i++)
  75     {
  76         const real *xi = x[index[i]];
  77         for (j = i+1; (j < nind); j++)
  78         {
  79             pbc_dx(&pbc, xi, x[index[j]], dx);
  80             temp2   = norm2(dx);
  81             d[i][j] = std::sqrt(temp2);
  82         }
  83     }
  84 }
  85
  86 static void calc_dist_tot(int nind, int index[], rvec x[],
  87                           int ePBC, matrix box,
  88                           real **d, real **dtot, real **dtot2,
  89                           gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
  90 {
  91     int      i, j;
  92     real    *xi;
  93     real     temp, temp2, temp1_3;
  94     rvec     dx;
  95     t_pbc    pbc;
  96
  97     set_pbc(&pbc, ePBC, box);
  98     for (i = 0; (i < nind-1); i++)
  99     {
 100         xi = x[index[i]];
 101         for (j = i+1; (j < nind); j++)
 102         {
 103             pbc_dx(&pbc, xi, x[index[j]], dx);
 104             temp2        = norm2(dx);
 105             temp         = std::sqrt(temp2);
 106             d[i][j]      = temp;
 107             dtot[i][j]  += temp;
 108             dtot2[i][j] += temp2;
 109             if (bNMR)
 110             {
 111                 temp1_3        = 1.0/(temp*temp2);
 112                 dtot1_3[i][j] += temp1_3;
 113                 dtot1_6[i][j] += temp1_3*temp1_3;
 114             }
 115         }
 116     }
 117 }
 118
 119 static void calc_nmr(int nind, int nframes, real **dtot1_3, real **dtot1_6,
 120                      real *max1_3, real *max1_6)
 121 {
 122     int     i, j;
 123     real    temp1_3, temp1_6;
 124
 125     for (i = 0; (i < nind-1); i++)
 126     {
 127         for (j = i+1; (j < nind); j++)
 128         {
 129             temp1_3 = gmx::invcbrt(dtot1_3[i][j]/nframes);
 130             temp1_6 = gmx::invsixthroot(dtot1_6[i][j]/nframes);
 131             if (temp1_3 > *max1_3)
 132             {
 133                 *max1_3 = temp1_3;
 134             }
 135             if (temp1_6 > *max1_6)
 136             {
 137                 *max1_6 = temp1_6;
 138             }
 139             dtot1_3[i][j] = temp1_3;
 140             dtot1_6[i][j] = temp1_6;
 141             dtot1_3[j][i] = temp1_3;
 142             dtot1_6[j][i] = temp1_6;
 143         }
 144     }
 145 }
 146
 147 static char Hnum[] = "123";
 148
 149 typedef struct {
 150     int  nr;
 151     real r_3;
 152     real r_6;
 153     real i_3;
 154     real i_6;
 155 } t_noe;
 156
 157 typedef struct {
 158     int   anr;
 159     int   ianr;
 160     int   rnr;
 161     char *aname;
 162     char *rname;
 163 } t_noe_gr;
 164
 165 typedef struct {
 166     bool  set;
 167     int   rnr;
 168     char *nname;
 169     char *rname;
 170     char *aname;
 171 } t_equiv;
 172
 173 static int read_equiv(const char *eq_fn, t_equiv ***equivptr)
 174 {
 175     FILE     *fp;
 176     char      line[STRLEN], resname[10], atomname[10], *lp;
 177     int       neq, na, n, resnr;
 178     t_equiv **equiv;
 179
 180     fp    = gmx_ffopen(eq_fn, "r");
 181     neq   = 0;
 182     equiv = nullptr;
 183     while (get_a_line(fp, line, STRLEN))
 184     {
 185         lp = line;
 186         /* this is not efficient, but I'm lazy */
 187         srenew(equiv, neq+1);
 188         equiv[neq] = nullptr;
 189         na         = 0;
 190         if (sscanf(lp, "%s %n", atomname, &n) == 1)
 191         {
 192             lp += n;
 193             snew(equiv[neq], 1);
 194             equiv[neq][0].nname = gmx_strdup(atomname);
 195             while (sscanf(lp, "%d %s %s %n", &resnr, resname, atomname, &n) == 3)
 196             {
 197                 /* this is not efficient, but I'm lazy (again) */
 198                 srenew(equiv[neq], na+1);
 199                 equiv[neq][na].set   = true;
 200                 equiv[neq][na].rnr   = resnr-1;
 201                 equiv[neq][na].rname = gmx_strdup(resname);
 202                 equiv[neq][na].aname = gmx_strdup(atomname);
 203                 if (na > 0)
 204                 {
 205                     equiv[neq][na].nname = nullptr;
 206                 }
 207                 na++;
 208                 lp += n;
 209             }
 210         }
 211         /* make empty element as flag for end of array */
 212         srenew(equiv[neq], na+1);
 213         equiv[neq][na].set   = false;
 214         equiv[neq][na].rnr   = 0;
 215         equiv[neq][na].rname = nullptr;
 216         equiv[neq][na].aname = nullptr;
 217
 218         /* next */
 219         neq++;
 220     }
 221     gmx_ffclose(fp);
 222
 223     *equivptr = equiv;
 224
 225     return neq;
 226 }
 227
 228 static void dump_equiv(FILE *out, int neq, t_equiv **equiv)
 229 {
 230     int i, j;
 231
 232     fprintf(out, "Dumping equivalent list\n");
 233     for (i = 0; i < neq; i++)
 234     {
 235         fprintf(out, "%s", equiv[i][0].nname);
 236         for (j = 0; equiv[i][j].set; j++)
 237         {
 238             fprintf(out, " %d %s %s",
 239                     equiv[i][j].rnr, equiv[i][j].rname, equiv[i][j].aname);
 240         }
 241         fprintf(out, "\n");
 242     }
 243 }
 244
 245 static gmx_bool is_equiv(int neq, t_equiv **equiv, char **nname,
 246                          int rnr1, char *rname1, char *aname1,
 247                          int rnr2, char *rname2, char *aname2)
 248 {
 249     int      i, j;
 250     gmx_bool bFound;
 251
 252     bFound = FALSE;
 253     /* we can terminate each loop when bFound is true! */
 254     for (i = 0; i < neq && !bFound; i++)
 255     {
 256         /* find first atom */
 257         for (j = 0; equiv[i][j].set && !bFound; j++)
 258         {
 259             bFound = ( equiv[i][j].rnr == rnr1 &&
 260                        std::strcmp(equiv[i][j].rname, rname1) == 0 &&
 261                        std::strcmp(equiv[i][j].aname, aname1) == 0 );
 262         }
 263         if (bFound)
 264         {
 265             /* find second atom */
 266             bFound = FALSE;
 267             for (j = 0; equiv[i][j].set && !bFound; j++)
 268             {
 269                 bFound = ( equiv[i][j].rnr == rnr2 &&
 270                            std::strcmp(equiv[i][j].rname, rname2) == 0 &&
 271                            std::strcmp(equiv[i][j].aname, aname2) == 0 );
 272             }
 273         }
 274     }
 275     if (bFound)
 276     {
 277         *nname = gmx_strdup(equiv[i-1][0].nname);
 278     }
 279
 280     return bFound;
 281 }
 282
 283 static int analyze_noe_equivalent(const char *eq_fn,
 284                                   const t_atoms *atoms, int isize, int *index,
 285                                   gmx_bool bSumH,
 286                                   int *noe_index, t_noe_gr *noe_gr)
 287 {
 288     int       i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
 289     char     *anmi, *anmj, **nnm;
 290     gmx_bool  bMatch, bEquiv;
 291     t_equiv **equiv;
 292
 293     snew(nnm, isize);
 294     if (bSumH)
 295     {
 296         if (eq_fn)
 297         {
 298             neq = read_equiv(eq_fn, &equiv);
 299             if (debug)
 300             {
 301                 dump_equiv(debug, neq, equiv);
 302             }
 303         }
 304         else
 305         {
 306             neq   = 0;
 307             equiv = nullptr;
 308         }
 309
 310         groupnr = 0;
 311         for (i = 0; i < isize; i++)
 312         {
 313             if (equiv && i < isize-1)
 314             {
 315                 /* check explicit list of equivalent atoms */
 316                 do
 317                 {
 318                     j      = i+1;
 319                     rnri   = atoms->atom[index[i]].resind;
 320                     rnrj   = atoms->atom[index[j]].resind;
 321                     bEquiv =
 322                         is_equiv(neq, equiv, &nnm[i],
 323                                  rnri, *atoms->resinfo[rnri].name, *atoms->atomname[index[i]],
 324                                  rnrj, *atoms->resinfo[rnrj].name, *atoms->atomname[index[j]]);
 325                     if (nnm[i] && bEquiv)
 326                     {
 327                         nnm[j] = gmx_strdup(nnm[i]);
 328                     }
 329                     if (bEquiv)
 330                     {
 331                         /* set index for matching atom */
 332                         noe_index[i] = groupnr;
 333                         /* skip matching atom */
 334                         i = j;
 335                     }
 336                 }
 337                 while (bEquiv && i < isize-1);
 338             }
 339             else
 340             {
 341                 bEquiv = FALSE;
 342             }
 343             if (!bEquiv)
 344             {
 345                 /* look for triplets of consecutive atoms with name XX?,
 346                    X are any number of letters or digits and ? goes from 1 to 3
 347                    This is supposed to cover all CH3 groups and the like */
 348                 anmi   = *atoms->atomname[index[i]];
 349                 anmil  = std::strlen(anmi);
 350                 bMatch = i <= isize-3 && anmi[anmil-1] == '1';
 351                 if (bMatch)
 352                 {
 353                     for (j = 1; j < 3; j++)
 354                     {
 355                         anmj   = *atoms->atomname[index[i+j]];
 356                         anmjl  = std::strlen(anmj);
 357                         bMatch = bMatch && ( anmil == anmjl && anmj[anmjl-1] == Hnum[j] &&
 358                                              std::strncmp(anmi, anmj, anmil-1) == 0 );
 359                     }
 360                 }
 361                 /* set index for this atom */
 362                 noe_index[i] = groupnr;
 363                 if (bMatch)
 364                 {
 365                     /* set index for next two matching atoms */
 366                     for (j = 1; j < 3; j++)
 367                     {
 368                         noe_index[i+j] = groupnr;
 369                     }
 370                     /* skip matching atoms */
 371                     i += 2;
 372                 }
 373             }
 374             groupnr++;
 375         }
 376     }
 377     else
 378     {
 379         /* make index without looking for equivalent atoms */
 380         for (i = 0; i < isize; i++)
 381         {
 382             noe_index[i] = i;
 383         }
 384         groupnr = isize;
 385     }
 386     noe_index[isize] = groupnr;
 387
 388     if (debug)
 389     {
 390         /* dump new names */
 391         for (i = 0; i < isize; i++)
 392         {
 393             rnri = atoms->atom[index[i]].resind;
 394             fprintf(debug, "%s %s %d -> %s\n", *atoms->atomname[index[i]],
 395                     *atoms->resinfo[rnri].name, rnri, nnm[i] ? nnm[i] : "");
 396         }
 397     }
 398
 399     for (i = 0; i < isize; i++)
 400     {
 401         gi = noe_index[i];
 402         if (!noe_gr[gi].aname)
 403         {
 404             noe_gr[gi].ianr = i;
 405             noe_gr[gi].anr  = index[i];
 406             if (nnm[i])
 407             {
 408                 noe_gr[gi].aname = gmx_strdup(nnm[i]);
 409             }
 410             else
 411             {
 412                 noe_gr[gi].aname = gmx_strdup(*atoms->atomname[index[i]]);
 413                 if (noe_index[i] == noe_index[i+1])
 414                 {
 415                     noe_gr[gi].aname[std::strlen(noe_gr[gi].aname)-1] = '*';
 416                 }
 417             }
 418             noe_gr[gi].rnr   = atoms->atom[index[i]].resind;
 419             noe_gr[gi].rname = gmx_strdup(*atoms->resinfo[noe_gr[gi].rnr].name);
 420             /* dump group definitions */
 421             if (debug)
 422             {
 423                 fprintf(debug, "%d %d %d %d %s %s %d\n", i, gi,
 424                         noe_gr[gi].ianr, noe_gr[gi].anr, noe_gr[gi].aname,
 425                         noe_gr[gi].rname, noe_gr[gi].rnr);
 426             }
 427         }
 428     }
 429     for (i = 0; i < isize; i++)
 430     {
 431         sfree(nnm[i]);
 432     }
 433     sfree(nnm);
 434
 435     return groupnr;
 436 }
 437
 438 /* #define NSCALE 3 */
 439 /*  static char *noe_scale[NSCALE+1] = { */
 440 /*    "strong", "medium", "weak", "none" */
 441 /*  }; */
 442 #define NSCALE 6
 443
 444 static char *noe2scale(real r3, real r6, real rmax)
 445 {
 446     int         i, s3, s6;
 447     static char buf[NSCALE+1];
 448
 449     /* r goes from 0 to rmax
 450        NSCALE*r/rmax goes from 0 to NSCALE
 451        NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
 452     s3 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE*r3/rmax));
 453     s6 = NSCALE - std::min(NSCALE, static_cast<int>(NSCALE*r6/rmax));
 454
 455     for (i = 0; i < s3; i++)
 456     {
 457         buf[i] = '=';
 458     }
 459     for (; i < s6; i++)
 460     {
 461         buf[i] = '-';
 462     }
 463     buf[i] = '\0';
 464
 465     return buf;
 466 }
 467
 468 static void calc_noe(int isize, int *noe_index,
 469                      real **dtot1_3, real **dtot1_6, int gnr, t_noe **noe)
 470 {
 471     int i, j, gi, gj;
 472
 473     /* make half matrix of noe-group distances from atom distances */
 474     for (i = 0; i < isize; i++)
 475     {
 476         gi = noe_index[i];
 477         for (j = i; j < isize; j++)
 478         {
 479             gj = noe_index[j];
 480             noe[gi][gj].nr++;
 481             noe[gi][gj].i_3 += 1.0/gmx::power3(dtot1_3[i][j]);
 482             noe[gi][gj].i_6 += 1.0/gmx::power6(dtot1_6[i][j]);
 483         }
 484     }
 485
 486     /* make averages */
 487     for (i = 0; i < gnr; i++)
 488     {
 489         for (j = i+1; j < gnr; j++)
 490         {
 491             noe[i][j].r_3 = gmx::invcbrt(noe[i][j].i_3/noe[i][j].nr);
 492             noe[i][j].r_6 = gmx::invsixthroot(noe[i][j].i_6/noe[i][j].nr);
 493             noe[j][i]     = noe[i][j];
 494         }
 495     }
 496 }
 497
 498 static void write_noe(FILE *fp, int gnr, t_noe **noe, t_noe_gr *noe_gr, real rmax)
 499 {
 500     int      i, j;
 501     real     r3, r6, min3, min6;
 502     char     buf[10], b3[10], b6[10];
 503     t_noe_gr gri, grj;
 504
 505     min3 = min6 = 1e6;
 506     fprintf(fp,
 507             ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
 508             "ianr", "anr", "anm", "rnm", "rnr", "ianr", "anr", "anm", "rnm", "rnr",
 509             "1/r^3", "1/r^6", "intnsty", "Dr", "Da", "scale");
 510     for (i = 0; i < gnr; i++)
 511     {
 512         gri = noe_gr[i];
 513         for (j = i+1; j < gnr; j++)
 514         {
 515             grj  = noe_gr[j];
 516             r3   = noe[i][j].r_3;
 517             r6   = noe[i][j].r_6;
 518             min3 = std::min(r3, min3);
 519             min6 = std::min(r6, min6);
 520             if (r3 < rmax || r6 < rmax)
 521             {
 522                 if (grj.rnr == gri.rnr)
 523                 {
 524                     sprintf(buf, "%2d", grj.anr-gri.anr);
 525                 }
 526                 else
 527                 {
 528                     buf[0] = '\0';
 529                 }
 530                 if (r3 < rmax)
 531                 {
 532                     sprintf(b3, "%-5.3f", r3);
 533                 }
 534                 else
 535                 {
 536                     std::strcpy(b3, "-");
 537                 }
 538                 if (r6 < rmax)
 539                 {
 540                     sprintf(b6, "%-5.3f", r6);
 541                 }
 542                 else
 543                 {
 544                     std::strcpy(b6, "-");
 545                 }
 546                 fprintf(fp,
 547                         "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
 548                         gri.ianr+1, gri.anr+1, gri.aname, gri.rname, gri.rnr+1,
 549                         grj.ianr+1, grj.anr+1, grj.aname, grj.rname, grj.rnr+1,
 550                         b3, b6, static_cast<int>(noe[i][j].i_6+0.5), grj.rnr-gri.rnr, buf,
 551                         noe2scale(r3, r6, rmax));
 552             }
 553         }
 554     }
 555 #define MINI ((i == 3) ? min3 : min6)
 556     for (i = 3; i <= 6; i += 3)
 557     {
 558         if (MINI > rmax)
 559         {
 560             fprintf(stdout, "NOTE: no 1/r^%d averaged distances found below %g, "
 561                     "smallest was %g\n", i, rmax, MINI );
 562         }
 563         else
 564         {
 565             fprintf(stdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI );
 566         }
 567     }
 568 #undef MINI
 569 }
 570
 571 static void calc_rms(int nind, int nframes,
 572                      real **dtot, real **dtot2,
 573                      real **rmsmat,  real *rmsmax,
 574                      real **rmscmat, real *rmscmax,
 575                      real **meanmat, real *meanmax)
 576 {
 577     int     i, j;
 578     real    mean, mean2, rms, rmsc;
 579 /* N.B. dtot and dtot2 contain the total distance and the total squared
 580  * distance respectively, BUT they return RMS and the scaled RMS resp.
 581  */
 582
 583     *rmsmax  = -1000;
 584     *rmscmax = -1000;
 585     *meanmax = -1000;
 586
 587     for (i = 0; (i < nind-1); i++)
 588     {
 589         for (j = i+1; (j < nind); j++)
 590         {
 591             mean  = dtot[i][j]/nframes;
 592             mean2 = dtot2[i][j]/nframes;
 593             rms   = std::sqrt(std::max(static_cast<real>(0.0), mean2-mean*mean));
 594             rmsc  = rms/mean;
 595             if (mean > *meanmax)
 596             {
 597                 *meanmax = mean;
 598             }
 599             if (rms  > *rmsmax)
 600             {
 601                 *rmsmax = rms;
 602             }
 603             if (rmsc > *rmscmax)
 604             {
 605                 *rmscmax = rmsc;
 606             }
 607             meanmat[i][j] = meanmat[j][i] = mean;
 608             rmsmat[i][j]  = rmsmat[j][i] = rms;
 609             rmscmat[i][j] = rmscmat[j][i] = rmsc;
 610         }
 611     }
 612 }
 613
 614 real rms_diff(int natom, real **d, real **d_r)
 615 {
 616     int  i, j;
 617     real r, r2;
 618
 619     r2 = 0.0;
 620     for (i = 0; (i < natom-1); i++)
 621     {
 622         for (j = i+1; (j < natom); j++)
 623         {
 624             r   = d[i][j]-d_r[i][j];
 625             r2 += r*r;
 626         }
 627     }
 628     r2 /= (natom*(natom-1))/2;
 629
 630     return std::sqrt(r2);
 631 }
 632
 633 int gmx_rmsdist(int argc, char *argv[])
 634 {
 635     const char       *desc[] = {
 636         "[THISMODULE] computes the root mean square deviation of atom distances,",
 637         "which has the advantage that no fit is needed like in standard RMS",
 638         "deviation as computed by [gmx-rms].",
 639         "The reference structure is taken from the structure file.",
 640         "The RMSD at time t is calculated as the RMS",
 641         "of the differences in distance between atom-pairs in the reference",
 642         "structure and the structure at time t.[PAR]",
 643         "[THISMODULE] can also produce matrices of the rms distances, rms distances",
 644         "scaled with the mean distance and the mean distances and matrices with",
 645         "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
 646         "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
 647         "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
 648         "can be generated, by default averaging over equivalent hydrogens",
 649         "(all triplets of hydrogens named \\*[123]). Additionally a list of",
 650         "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
 651         "a set of equivalent atoms specified as residue number and name and",
 652         "atom name; e.g.:[PAR]",
 653         "[TT]HB* 3 SER  HB1 3 SER  HB2[tt][PAR]",
 654         "Residue and atom names must exactly match those in the structure",
 655         "file, including case. Specifying non-sequential atoms is undefined."
 656
 657     };
 658
 659     int               i, teller;
 660     real              t;
 661
 662     t_topology        top;
 663     int               ePBC;
 664     t_atoms          *atoms;
 665     matrix            box;
 666     rvec             *x;
 667     FILE             *fp;
 668
 669     t_trxstatus      *status;
 670     int               isize, gnr = 0;
 671     int              *index, *noe_index;
 672     char             *grpname;
 673     real            **d_r, **d, **dtot, **dtot2, **mean, **rms, **rmsc, *resnr;
 674     real            **dtot1_3 = nullptr, **dtot1_6 = nullptr;
 675     real              rmsnow, meanmax, rmsmax, rmscmax;
 676     real              max1_3, max1_6;
 677     t_noe_gr         *noe_gr = nullptr;
 678     t_noe           **noe    = nullptr;
 679     t_rgb             rlo, rhi;
 680     gmx_bool          bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
 681
 682     static int        nlevels  = 40;
 683     static real       scalemax = -1.0;
 684     static gmx_bool   bSumH    = TRUE;
 685     static gmx_bool   bPBC     = TRUE;
 686     gmx_output_env_t *oenv;
 687
 688     t_pargs           pa[] = {
 689         { "-nlevels",   FALSE, etINT,  {&nlevels},
 690           "Discretize RMS in this number of levels" },
 691         { "-max",   FALSE, etREAL, {&scalemax},
 692           "Maximum level in matrices" },
 693         { "-sumh",  FALSE, etBOOL, {&bSumH},
 694           "Average distance over equivalent hydrogens" },
 695         { "-pbc",   FALSE, etBOOL, {&bPBC},
 696           "Use periodic boundary conditions when computing distances" }
 697     };
 698     t_filenm          fnm[] = {
 699         { efTRX, "-f",   nullptr,       ffREAD },
 700         { efTPS, nullptr,   nullptr,       ffREAD },
 701         { efNDX, nullptr,   nullptr,       ffOPTRD },
 702         { efDAT, "-equiv", "equiv",   ffOPTRD },
 703         { efXVG, nullptr,   "distrmsd", ffWRITE },
 704         { efXPM, "-rms", "rmsdist",  ffOPTWR },
 705         { efXPM, "-scl", "rmsscale", ffOPTWR },
 706         { efXPM, "-mean", "rmsmean",  ffOPTWR },
 707         { efXPM, "-nmr3", "nmr3",     ffOPTWR },
 708         { efXPM, "-nmr6", "nmr6",     ffOPTWR },
 709         { efDAT, "-noe", "noe",     ffOPTWR },
 710     };
 711 #define NFILE asize(fnm)
 712
 713     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 714                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 715     {
 716         return 0;
 717     }
 718
 719     bRMS   = opt2bSet("-rms", NFILE, fnm);
 720     bScale = opt2bSet("-scl", NFILE, fnm);
 721     bMean  = opt2bSet("-mean", NFILE, fnm);
 722     bNOE   = opt2bSet("-noe", NFILE, fnm);
 723     bNMR3  = opt2bSet("-nmr3", NFILE, fnm);
 724     bNMR6  = opt2bSet("-nmr6", NFILE, fnm);
 725     bNMR   = bNMR3 || bNMR6 || bNOE;
 726
 727     max1_3 = 0;
 728     max1_6 = 0;
 729
 730     /* check input */
 731     if (bNOE && scalemax < 0)
 732     {
 733         scalemax = 0.6;
 734         fprintf(stderr, "WARNING: using -noe without -max "
 735                 "makes no sense, setting -max to %g\n\n", scalemax);
 736     }
 737
 738     /* get topology and index */
 739     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box, FALSE);
 740
 741     if (!bPBC)
 742     {
 743         ePBC = epbcNONE;
 744     }
 745     atoms = &(top.atoms);
 746
 747     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 748
 749     /* initialize arrays */
 750     snew(d, isize);
 751     snew(dtot, isize);
 752     snew(dtot2, isize);
 753     if (bNMR)
 754     {
 755         snew(dtot1_3, isize);
 756         snew(dtot1_6, isize);
 757     }
 758     snew(mean, isize);
 759     snew(rms, isize);
 760     snew(rmsc, isize);
 761     snew(d_r, isize);
 762     snew(resnr, isize);
 763     for (i = 0; (i < isize); i++)
 764     {
 765         snew(d[i], isize);
 766         snew(dtot[i], isize);
 767         snew(dtot2[i], isize);
 768         if (bNMR)
 769         {
 770             snew(dtot1_3[i], isize);
 771             snew(dtot1_6[i], isize);
 772         }
 773         snew(mean[i], isize);
 774         snew(rms[i], isize);
 775         snew(rmsc[i], isize);
 776         snew(d_r[i], isize);
 777         resnr[i] = i+1;
 778     }
 779
 780     /*set box type*/
 781     calc_dist(isize, index, x, ePBC, box, d_r);
 782     sfree(x);
 783
 784     /*open output files*/
 785     fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "RMS Deviation", "Time (ps)", "RMSD (nm)",
 786                   oenv);
 787     if (output_env_get_print_xvgr_codes(oenv))
 788     {
 789         fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
 790     }
 791
 792     /*do a first step*/
 793     read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 794     teller = 0;
 795
 796     do
 797     {
 798         calc_dist_tot(isize, index, x, ePBC, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
 799
 800         rmsnow = rms_diff(isize, d, d_r);
 801         fprintf(fp, "%g  %g\n", t, rmsnow);
 802         teller++;
 803     }
 804     while (read_next_x(oenv, status, &t, x, box));
 805     fprintf(stderr, "\n");
 806
 807     xvgrclose(fp);
 808
 809     close_trj(status);
 810
 811     calc_rms(isize, teller, dtot, dtot2, rms, &rmsmax, rmsc, &rmscmax, mean, &meanmax);
 812     fprintf(stderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
 813
 814     if (bNMR)
 815     {
 816         calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
 817     }
 818
 819     if (scalemax > -1.0)
 820     {
 821         rmsmax  = scalemax;
 822         rmscmax = scalemax;
 823         meanmax = scalemax;
 824         max1_3  = scalemax;
 825         max1_6  = scalemax;
 826     }
 827
 828     if (bNOE)
 829     {
 830         /* make list of noe atom groups */
 831         snew(noe_index, isize+1);
 832         snew(noe_gr, isize);
 833         gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE, fnm),
 834                                      atoms, isize, index, bSumH, noe_index, noe_gr);
 835         fprintf(stdout, "Found %d non-equivalent atom-groups in %d atoms\n",
 836                 gnr, isize);
 837         /* make half matrix of of noe-group distances from atom distances */
 838         snew(noe, gnr);
 839         for (i = 0; i < gnr; i++)
 840         {
 841             snew(noe[i], gnr);
 842         }
 843         calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
 844     }
 845
 846     rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
 847     rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;
 848
 849     if (bRMS)
 850     {
 851         write_xpm(opt2FILE("-rms", NFILE, fnm, "w"), 0,
 852                   "RMS of distance", "RMS (nm)", "Atom Index", "Atom Index",
 853                   isize, isize, resnr, resnr, rms, 0.0, rmsmax, rlo, rhi, &nlevels);
 854     }
 855
 856     if (bScale)
 857     {
 858         write_xpm(opt2FILE("-scl", NFILE, fnm, "w"), 0,
 859                   "Relative RMS", "RMS", "Atom Index", "Atom Index",
 860                   isize, isize, resnr, resnr, rmsc, 0.0, rmscmax, rlo, rhi, &nlevels);
 861     }
 862
 863     if (bMean)
 864     {
 865         write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0,
 866                   "Mean Distance", "Distance (nm)", "Atom Index", "Atom Index",
 867                   isize, isize, resnr, resnr, mean, 0.0, meanmax, rlo, rhi, &nlevels);
 868     }
 869
 870     if (bNMR3)
 871     {
 872         write_xpm(opt2FILE("-nmr3", NFILE, fnm, "w"), 0, "1/r^3 averaged distances",
 873                   "Distance (nm)", "Atom Index", "Atom Index",
 874                   isize, isize, resnr, resnr, dtot1_3, 0.0, max1_3, rlo, rhi, &nlevels);
 875     }
 876     if (bNMR6)
 877     {
 878         write_xpm(opt2FILE("-nmr6", NFILE, fnm, "w"), 0, "1/r^6 averaged distances",
 879                   "Distance (nm)", "Atom Index", "Atom Index",
 880                   isize, isize, resnr, resnr, dtot1_6, 0.0, max1_6, rlo, rhi, &nlevels);
 881     }
 882
 883     if (bNOE)
 884     {
 885         write_noe(opt2FILE("-noe", NFILE, fnm, "w"), gnr, noe, noe_gr, scalemax);
 886     }
 887
 888     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);
 889
 890     return 0;
 891 }