src/gromacs/gmxana/gmx_rmsdist.c

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