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[], int ePBC, matrix box, 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,
  87                           const int index[],
  88                           rvec      x[],
  89                           int       ePBC,
  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, ePBC, 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, 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)
 333                     {
 334                         nnm[j] = gmx_strdup(nnm[i]);
 335                     }
 336                     if (bEquiv)
 337                     {
 338                         /* set index for matching atom */
 339                         noe_index[i] = groupnr;
 340                         /* skip matching atom */
 341                         i = j;
 342                     }
 343                 } while (bEquiv && i < isize - 1);
 344             }
 345             else
 346             {
 347                 bEquiv = FALSE;
 348             }
 349             if (!bEquiv)
 350             {
 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';
 357                 if (bMatch)
 358                 {
 359                     for (j = 1; j < 3; j++)
 360                     {
 361                         anmj   = *atoms->atomname[index[i + j]];
 362                         anmjl  = std::strlen(anmj);
 363                         bMatch = bMatch
 364                                  && (anmil == anmjl && anmj[anmjl - 1] == Hnum[j]
 365                                      && std::strncmp(anmi, anmj, anmil - 1) == 0);
 366                     }
 367                 }
 368                 /* set index for this atom */
 369                 noe_index[i] = groupnr;
 370                 if (bMatch)
 371                 {
 372                     /* set index for next two matching atoms */
 373                     for (j = 1; j < 3; j++)
 374                     {
 375                         noe_index[i + j] = groupnr;
 376                     }
 377                     /* skip matching atoms */
 378                     i += 2;
 379                 }
 380             }
 381             groupnr++;
 382         }
 383     }
 384     else
 385     {
 386         /* make index without looking for equivalent atoms */
 387         for (i = 0; i < isize; i++)
 388         {
 389             noe_index[i] = i;
 390         }
 391         groupnr = isize;
 392     }
 393     noe_index[isize] = groupnr;
 394
 395     if (debug)
 396     {
 397         /* dump new names */
 398         for (i = 0; i < isize; i++)
 399         {
 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] : "");
 403         }
 404     }
 405
 406     for (i = 0; i < isize; i++)
 407     {
 408         gi = noe_index[i];
 409         if (!noe_gr[gi].aname)
 410         {
 411             noe_gr[gi].ianr = i;
 412             noe_gr[gi].anr  = index[i];
 413             if (nnm[i])
 414             {
 415                 noe_gr[gi].aname = gmx_strdup(nnm[i]);
 416             }
 417             else
 418             {
 419                 noe_gr[gi].aname = gmx_strdup(*atoms->atomname[index[i]]);
 420                 if (noe_index[i] == noe_index[i + 1])
 421                 {
 422                     noe_gr[gi].aname[std::strlen(noe_gr[gi].aname) - 1] = '*';
 423                 }
 424             }
 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 */
 428             if (debug)
 429             {
 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);
 432             }
 433         }
 434     }
 435     for (i = 0; i < isize; i++)
 436     {
 437         sfree(nnm[i]);
 438     }
 439     sfree(nnm);
 440
 441     return groupnr;
 442 }
 443
 444 /* #define NSCALE 3 */
 445 /*  static char *noe_scale[NSCALE+1] = { */
 446 /*    "strong", "medium", "weak", "none" */
 447 /*  }; */
 448 #define NSCALE 6
 449
 450 static char* noe2scale(real r3, real r6, real rmax)
 451 {
 452     int         i, s3, s6;
 453     static char buf[NSCALE + 1];
 454
 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));
 460
 461     for (i = 0; i < s3; i++)
 462     {
 463         buf[i] = '=';
 464     }
 465     for (; i < s6; i++)
 466     {
 467         buf[i] = '-';
 468     }
 469     buf[i] = '\0';
 470
 471     return buf;
 472 }
 473
 474 static void calc_noe(int isize, const int* noe_index, real** dtot1_3, real** dtot1_6, int gnr, t_noe** noe)
 475 {
 476     int i, j, gi, gj;
 477
 478     /* make half matrix of noe-group distances from atom distances */
 479     for (i = 0; i < isize; i++)
 480     {
 481         gi = noe_index[i];
 482         for (j = i; j < isize; j++)
 483         {
 484             gj = noe_index[j];
 485             noe[gi][gj].nr++;
 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]);
 488         }
 489     }
 490
 491     /* make averages */
 492     for (i = 0; i < gnr; i++)
 493     {
 494         for (j = i + 1; j < gnr; j++)
 495         {
 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];
 499         }
 500     }
 501 }
 502
 503 static void write_noe(FILE* fp, int gnr, t_noe** noe, t_noe_gr* noe_gr, real rmax)
 504 {
 505     int      i, j;
 506     real     r3, r6, min3, min6;
 507     char     buf[10], b3[10], b6[10];
 508     t_noe_gr gri, grj;
 509
 510     min3 = min6 = 1e6;
 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++)
 515     {
 516         gri = noe_gr[i];
 517         for (j = i + 1; j < gnr; j++)
 518         {
 519             grj  = noe_gr[j];
 520             r3   = noe[i][j].r_3;
 521             r6   = noe[i][j].r_6;
 522             min3 = std::min(r3, min3);
 523             min6 = std::min(r6, min6);
 524             if (r3 < rmax || r6 < rmax)
 525             {
 526                 if (grj.rnr == gri.rnr)
 527                 {
 528                     sprintf(buf, "%2d", grj.anr - gri.anr);
 529                 }
 530                 else
 531                 {
 532                     buf[0] = '\0';
 533                 }
 534                 if (r3 < rmax)
 535                 {
 536                     sprintf(b3, "%-5.3f", r3);
 537                 }
 538                 else
 539                 {
 540                     std::strcpy(b3, "-");
 541                 }
 542                 if (r6 < rmax)
 543                 {
 544                     sprintf(b6, "%-5.3f", r6);
 545                 }
 546                 else
 547                 {
 548                     std::strcpy(b6, "-");
 549                 }
 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));
 554             }
 555         }
 556     }
 557 #define MINI ((i == 3) ? min3 : min6)
 558     for (i = 3; i <= 6; i += 3)
 559     {
 560         if (MINI > rmax)
 561         {
 562             fprintf(stdout,
 563                     "NOTE: no 1/r^%d averaged distances found below %g, "
 564                     "smallest was %g\n",
 565                     i, rmax, MINI);
 566         }
 567         else
 568         {
 569             fprintf(stdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI);
 570         }
 571     }
 572 #undef MINI
 573 }
 574
 575 static void calc_rms(int    nind,
 576                      int    nframes,
 577                      real** dtot,
 578                      real** dtot2,
 579                      real** rmsmat,
 580                      real*  rmsmax,
 581                      real** rmscmat,
 582                      real*  rmscmax,
 583                      real** meanmat,
 584                      real*  meanmax)
 585 {
 586     int  i, j;
 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.
 590      */
 591
 592     *rmsmax  = -1000;
 593     *rmscmax = -1000;
 594     *meanmax = -1000;
 595
 596     for (i = 0; (i < nind - 1); i++)
 597     {
 598         for (j = i + 1; (j < nind); j++)
 599         {
 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));
 603             rmsc  = rms / mean;
 604             if (mean > *meanmax)
 605             {
 606                 *meanmax = mean;
 607             }
 608             if (rms > *rmsmax)
 609             {
 610                 *rmsmax = rms;
 611             }
 612             if (rmsc > *rmscmax)
 613             {
 614                 *rmscmax = rmsc;
 615             }
 616             meanmat[i][j] = meanmat[j][i] = mean;
 617             rmsmat[i][j] = rmsmat[j][i] = rms;
 618             rmscmat[i][j] = rmscmat[j][i] = rmsc;
 619         }
 620     }
 621 }
 622
 623 static real rms_diff(int natom, real** d, real** d_r)
 624 {
 625     int  i, j;
 626     real r, r2;
 627
 628     r2 = 0.0;
 629     for (i = 0; (i < natom - 1); i++)
 630     {
 631         for (j = i + 1; (j < natom); j++)
 632         {
 633             r = d[i][j] - d_r[i][j];
 634             r2 += r * r;
 635         }
 636     }
 637     r2 /= gmx::exactDiv(natom * (natom - 1), 2);
 638
 639     return std::sqrt(r2);
 640 }
 641
 642 int gmx_rmsdist(int argc, char* argv[])
 643 {
 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."
 665
 666     };
 667
 668     int  i, teller;
 669     real t;
 670
 671     t_topology top;
 672     int        ePBC;
 673     t_atoms*   atoms;
 674     matrix     box;
 675     rvec*      x;
 676     FILE*      fp;
 677
 678     t_trxstatus* status;
 679     int          isize, gnr = 0;
 680     int *        index, *noe_index;
 681     char*        grpname;
 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;
 685     real         max1_3, max1_6;
 686     t_noe_gr*    noe_gr = nullptr;
 687     t_noe**      noe    = nullptr;
 688     t_rgb        rlo, rhi;
 689     gmx_bool     bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
 690
 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;
 696
 697     t_pargs pa[] = {
 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" },
 701         { "-pbc",
 702           FALSE,
 703           etBOOL,
 704           { &bPBC },
 705           "Use periodic boundary conditions when computing distances" }
 706     };
 707     t_filenm fnm[] = {
 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 },
 714     };
 715 #define NFILE asize(fnm)
 716
 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))
 719     {
 720         return 0;
 721     }
 722
 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;
 730
 731     max1_3 = 0;
 732     max1_6 = 0;
 733
 734     /* check input */
 735     if (bNOE && scalemax < 0)
 736     {
 737         scalemax = 0.6;
 738         fprintf(stderr,
 739                 "WARNING: using -noe without -max "
 740                 "makes no sense, setting -max to %g\n\n",
 741                 scalemax);
 742     }
 743
 744     /* get topology and index */
 745     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box, FALSE);
 746
 747     if (!bPBC)
 748     {
 749         ePBC = epbcNONE;
 750     }
 751     atoms = &(top.atoms);
 752
 753     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 754
 755     /* initialize arrays */
 756     snew(d, isize);
 757     snew(dtot, isize);
 758     snew(dtot2, isize);
 759     if (bNMR)
 760     {
 761         snew(dtot1_3, isize);
 762         snew(dtot1_6, isize);
 763     }
 764     snew(mean, isize);
 765     snew(rms, isize);
 766     snew(rmsc, isize);
 767     snew(d_r, isize);
 768     snew(resnr, isize);
 769     for (i = 0; (i < isize); i++)
 770     {
 771         snew(d[i], isize);
 772         snew(dtot[i], isize);
 773         snew(dtot2[i], isize);
 774         if (bNMR)
 775         {
 776             snew(dtot1_3[i], isize);
 777             snew(dtot1_6[i], isize);
 778         }
 779         snew(mean[i], isize);
 780         snew(rms[i], isize);
 781         snew(rmsc[i], isize);
 782         snew(d_r[i], isize);
 783         resnr[i] = i + 1;
 784     }
 785
 786     /*set box type*/
 787     calc_dist(isize, index, x, ePBC, box, d_r);
 788     sfree(x);
 789
 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))
 793     {
 794         fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
 795     }
 796
 797     /*do a first step*/
 798     read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 799     teller = 0;
 800
 801     do
 802     {
 803         calc_dist_tot(isize, index, x, ePBC, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
 804
 805         rmsnow = rms_diff(isize, d, d_r);
 806         fprintf(fp, "%g  %g\n", t, rmsnow);
 807         teller++;
 808     } while (read_next_x(oenv, status, &t, x, box));
 809     fprintf(stderr, "\n");
 810
 811     xvgrclose(fp);
 812
 813     close_trx(status);
 814
 815     calc_rms(isize, teller, dtot, dtot2, rms, &rmsmax, rmsc, &rmscmax, mean, &meanmax);
 816     fprintf(stderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
 817
 818     if (bNMR)
 819     {
 820         calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
 821     }
 822
 823     if (scalemax > -1.0)
 824     {
 825         rmsmax  = scalemax;
 826         rmscmax = scalemax;
 827         meanmax = scalemax;
 828         max1_3  = scalemax;
 829         max1_6  = scalemax;
 830     }
 831
 832     if (bNOE)
 833     {
 834         /* make list of noe atom groups */
 835         snew(noe_index, isize + 1);
 836         snew(noe_gr, isize);
 837         gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE, fnm), atoms, isize, index, bSumH,
 838                                      noe_index, noe_gr);
 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 */
 841         snew(noe, gnr);
 842         for (i = 0; i < gnr; i++)
 843         {
 844             snew(noe[i], gnr);
 845         }
 846         calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
 847     }
 848
 849     rlo.r = 1.0;
 850     rlo.g = 1.0;
 851     rlo.b = 1.0;
 852     rhi.r = 0.0;
 853     rhi.g = 0.0;
 854     rhi.b = 0.0;
 855
 856     if (bRMS)
 857     {
 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);
 860     }
 861
 862     if (bScale)
 863     {
 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);
 866     }
 867
 868     if (bMean)
 869     {
 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,
 872                   rhi, &nlevels);
 873     }
 874
 875     if (bNMR3)
 876     {
 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);
 880     }
 881     if (bNMR6)
 882     {
 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);
 886     }
 887
 888     if (bNOE)
 889     {
 890         write_noe(opt2FILE("-noe", NFILE, fnm, "w"), gnr, noe, noe_gr, scalemax);
 891     }
 892
 893     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), nullptr);
 894
 895     return 0;
 896 }