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