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