src/gromacs/gmxana/gmx_rmsdist.c

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