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