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