src/tools/gmx_mindist.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 <stdlib.h>
  41
  42 #include "sysstuff.h"
  43 #include "string.h"
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "xvgr.h"
  49 #include "pbc.h"
  50 #include "copyrite.h"
  51 #include "futil.h"
  52 #include "statutil.h"
  53 #include "index.h"
  54 #include "tpxio.h"
  55 #include "rmpbc.h"
  56 #include "xtcio.h"
  57 #include "gmx_ana.h"
  58
  59
  60 static void periodic_dist(matrix box,rvec x[],int n,atom_id index[],
  61                           real *rmin,real *rmax,int *min_ind)
  62 {
  63 #define NSHIFT 26
  64   int  sx,sy,sz,i,j,s;
  65   real sqr_box,r2min,r2max,r2;
  66   rvec shift[NSHIFT],d0,d;
  67
  68   sqr_box = sqr(min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ])));
  69
  70   s = 0;
  71   for(sz=-1; sz<=1; sz++)
  72     for(sy=-1; sy<=1; sy++)
  73       for(sx=-1; sx<=1; sx++)
  74         if (sx!=0 || sy!=0 || sz!=0) {
  75           for(i=0; i<DIM; i++)
  76             shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
  77           s++;
  78         }
  79
  80   r2min = sqr_box;
  81   r2max = 0;
  82
  83   for(i=0; i<n; i++)
  84     for(j=i+1; j<n; j++) {
  85       rvec_sub(x[index[i]],x[index[j]],d0);
  86       r2 = norm2(d0);
  87       if (r2 > r2max)
  88         r2max = r2;
  89       for(s=0; s<NSHIFT; s++) {
  90         rvec_add(d0,shift[s],d);
  91         r2 = norm2(d);
  92         if (r2 < r2min) {
  93           r2min = r2;
  94           min_ind[0] = i;
  95           min_ind[1] = j;
  96         }
  97       }
  98     }
  99
 100   *rmin = sqrt(r2min);
 101   *rmax = sqrt(r2max);
 102 }
 103
 104 static void periodic_mindist_plot(const char *trxfn,const char *outfn,
 105                                   t_topology *top,int ePBC,
 106                                   int n,atom_id index[],gmx_bool bSplit,
 107                                   const output_env_t oenv)
 108 {
 109   FILE   *out;
 110   const char *leg[5] = { "min per.","max int.","box1","box2","box3" };
 111   t_trxstatus *status;
 112   real   t;
 113   rvec   *x;
 114   matrix box;
 115   int    natoms,ind_min[2]={0,0},ind_mini=0,ind_minj=0;
 116   real   r,rmin,rmax,rmint,tmint;
 117   gmx_bool   bFirst;
 118   gmx_rmpbc_t  gpbc=NULL;
 119
 120   natoms=read_first_x(oenv,&status,trxfn,&t,&x,box);
 121
 122   check_index(NULL,n,index,NULL,natoms);
 123
 124   out = xvgropen(outfn,"Minimum distance to periodic image",
 125                  output_env_get_time_label(oenv),"Distance (nm)",oenv);
 126   if (output_env_get_print_xvgr_codes(oenv))
 127     fprintf(out,"@ subtitle \"and maximum internal distance\"\n");
 128   xvgr_legend(out,5,leg,oenv);
 129
 130   rmint = box[XX][XX];
 131   tmint = 0;
 132
 133   if (NULL != top)
 134     gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
 135
 136   bFirst=TRUE;
 137   do {
 138     if (NULL != top)
 139       gmx_rmpbc(gpbc,natoms,box,x);
 140
 141     periodic_dist(box,x,n,index,&rmin,&rmax,ind_min);
 142     if (rmin < rmint) {
 143       rmint = rmin;
 144       tmint = t;
 145       ind_mini = ind_min[0];
 146       ind_minj = ind_min[1];
 147     }
 148     if ( bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv))<1e-5 )
 149       fprintf(out, "&\n");
 150     fprintf(out,"\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
 151             output_env_conv_time(oenv,t),rmin,rmax,norm(box[0]),norm(box[1]),norm(box[2]));
 152     bFirst=FALSE;
 153   } while(read_next_x(oenv,status,&t,natoms,x,box));
 154
 155   if (NULL != top)
 156     gmx_rmpbc_done(gpbc);
 157
 158   ffclose(out);
 159
 160   fprintf(stdout,
 161           "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
 162           "between atoms %d and %d\n",
 163           rmint,output_env_conv_time(oenv,tmint),output_env_get_time_unit(oenv),
 164           index[ind_mini]+1,index[ind_minj]+1);
 165 }
 166
 167 static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
 168                       int nx1,int nx2, atom_id index1[], atom_id index2[],
 169                       gmx_bool bGroup,
 170                       real *rmin, real *rmax, int *nmin, int *nmax,
 171                       int *ixmin, int *jxmin, int *ixmax, int *jxmax)
 172 {
 173   int     i,j,i0=0,j1;
 174   int     ix,jx;
 175   atom_id *index3;
 176   rvec    dx;
 177   real    r2,rmin2,rmax2,rcut2;
 178   t_pbc   pbc;
 179   int     nmin_j,nmax_j;
 180
 181   *ixmin = -1;
 182   *jxmin = -1;
 183   *ixmax = -1;
 184   *jxmax = -1;
 185   *nmin = 0;
 186   *nmax = 0;
 187
 188   rcut2=sqr(rcut);
 189
 190   /* Must init pbc every step because of pressure coupling */
 191   if (bPBC)
 192     set_pbc(&pbc,ePBC,box);
 193   if (index2) {
 194     i0=0;
 195     j1=nx2;
 196     index3=index2;
 197   } else {
 198     j1=nx1;
 199     index3=index1;
 200   }
 201
 202   rmin2=1e12;
 203   rmax2=-1e12;
 204
 205   for(j=0; (j < j1); j++) {
 206     jx = index3[j];
 207     if (index2 == NULL) {
 208       i0 = j + 1;
 209     }
 210     nmin_j = 0;
 211     nmax_j = 0;
 212     for(i=i0; (i < nx1); i++) {
 213       ix = index1[i];
 214       if (ix != jx) {
 215         if (bPBC)
 216           pbc_dx(&pbc,x[ix],x[jx],dx);
 217         else
 218           rvec_sub(x[ix],x[jx],dx);
 219         r2=iprod(dx,dx);
 220         if (r2 < rmin2) {
 221           rmin2=r2;
 222           *ixmin=ix;
 223           *jxmin=jx;
 224         }
 225         if (r2 > rmax2) {
 226           rmax2=r2;
 227           *ixmax=ix;
 228           *jxmax=jx;
 229         }
 230         if (r2 <= rcut2) {
 231           nmin_j++;
 232         } else if (r2 > rcut2) {
 233           nmax_j++;
 234         }
 235       }
 236     }
 237     if (bGroup) {
 238       if (nmin_j > 0) {
 239         (*nmin)++;
 240       }
 241       if (nmax_j > 0) {
 242         (*nmax)++;
 243       }
 244     } else {
 245       *nmin += nmin_j;
 246       *nmax += nmax_j;
 247     }
 248   }
 249   *rmin = sqrt(rmin2);
 250   *rmax = sqrt(rmax2);
 251 }
 252
 253 void dist_plot(const char *fn,const char *afile,const char *dfile,
 254                const char *nfile,const char *rfile,const char *xfile,
 255                real rcut,gmx_bool bMat,t_atoms *atoms,
 256                int ng,atom_id *index[],int gnx[],char *grpn[],gmx_bool bSplit,
 257                gmx_bool bMin, int nres, atom_id *residue,gmx_bool bPBC,int ePBC,
 258                gmx_bool bGroup,gmx_bool bEachResEachTime, gmx_bool bPrintResName,
 259                const output_env_t oenv)
 260 {
 261   FILE         *atm,*dist,*num;
 262   t_trxstatus  *trxout;
 263   char         buf[256];
 264   char         **leg;
 265   real         t,dmin,dmax,**mindres=NULL,**maxdres=NULL;
 266   int          nmin,nmax;
 267   t_trxstatus  *status;
 268   int          i=-1,j,k,natoms;
 269   int          min1,min2,max1,max2,min1r,min2r,max1r,max2r;
 270   atom_id      oindex[2];
 271   rvec         *x0;
 272   matrix       box;
 273   t_trxframe   frout;
 274   gmx_bool         bFirst;
 275   FILE *respertime=NULL;
 276
 277   if ((natoms=read_first_x(oenv,&status,fn,&t,&x0,box))==0)
 278     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 279
 280   sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
 281   dist= xvgropen(dfile,buf,output_env_get_time_label(oenv),"Distance (nm)",oenv);
 282   sprintf(buf,"Number of Contacts %s %g nm",bMin ? "<" : ">",rcut);
 283   num = nfile ? xvgropen(nfile,buf,output_env_get_time_label(oenv),"Number",oenv) : NULL;
 284   atm = afile ? ffopen(afile,"w") : NULL;
 285   trxout = xfile ? open_trx(xfile,"w") : NULL;
 286
 287   if (bMat) {
 288     if (ng == 1) {
 289       snew(leg,1);
 290       sprintf(buf,"Internal in %s",grpn[0]);
 291       leg[0]=strdup(buf);
 292       xvgr_legend(dist,0,(const char**)leg,oenv);
 293       if (num) xvgr_legend(num,0,(const char**)leg,oenv);
 294     }
 295     else {
 296       snew(leg,(ng*(ng-1))/2);
 297       for(i=j=0; (i<ng-1); i++) {
 298         for(k=i+1; (k<ng); k++,j++) {
 299           sprintf(buf,"%s-%s",grpn[i],grpn[k]);
 300           leg[j]=strdup(buf);
 301         }
 302       }
 303       xvgr_legend(dist,j,(const char**)leg,oenv);
 304       if (num) xvgr_legend(num,j,(const char**)leg,oenv);
 305     }
 306   }
 307   else {
 308     snew(leg,ng-1);
 309     for(i=0; (i<ng-1); i++) {
 310       sprintf(buf,"%s-%s",grpn[0],grpn[i+1]);
 311       leg[i]=strdup(buf);
 312     }
 313     xvgr_legend(dist,ng-1,(const char**)leg,oenv);
 314     if (num) xvgr_legend(num,ng-1,(const char**)leg,oenv);
 315   }
 316
 317   if (bEachResEachTime)
 318   {
 319     sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
 320     respertime=xvgropen(rfile,buf,output_env_get_time_label(oenv),"Distance (nm)",oenv);
 321     xvgr_legend(respertime,ng-1,(const char**)leg,oenv);
 322         if (bPrintResName)
 323           fprintf(respertime,"# ");
 324             for (j=0; j<nres; j++)
 325                   fprintf(respertime,"%s%d ",*(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name),atoms->atom[index[0][residue[j]]].resind);
 326                 fprintf(respertime,"\n");
 327
 328   }
 329
 330   j=0;
 331   if (nres) {
 332     snew(mindres, ng-1);
 333     snew(maxdres, ng-1);
 334     for(i=1; i<ng; i++) {
 335       snew(mindres[i-1], nres);
 336       snew(maxdres[i-1], nres);
 337       for(j=0; j<nres; j++)
 338         mindres[i-1][j]=1e6;
 339       /* maxdres[*][*] is already 0 */
 340     }
 341   }
 342   bFirst=TRUE;
 343   do {
 344     if ( bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv))<1e-5 ) {
 345       fprintf(dist, "&\n");
 346       if (num) fprintf(num, "&\n");
 347       if (atm) fprintf(atm, "&\n");
 348     }
 349     fprintf(dist,"%12e",output_env_conv_time(oenv,t));
 350     if (num) fprintf(num,"%12e",output_env_conv_time(oenv,t));
 351
 352     if (bMat) {
 353       if (ng == 1) {
 354         calc_dist(rcut,bPBC,ePBC,box,x0,gnx[0],gnx[0],index[0],index[0],bGroup,
 355                   &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 356         fprintf(dist,"  %12e",bMin?dmin:dmax);
 357         if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 358       }
 359       else {
 360         for(i=0; (i<ng-1); i++) {
 361           for(k=i+1; (k<ng); k++) {
 362             calc_dist(rcut,bPBC,ePBC,box,x0,gnx[i],gnx[k],index[i],index[k],
 363                       bGroup,&dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 364             fprintf(dist,"  %12e",bMin?dmin:dmax);
 365             if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 366           }
 367         }
 368       }
 369     }
 370     else {
 371       for(i=1; (i<ng); i++) {
 372         calc_dist(rcut,bPBC,ePBC,box,x0,gnx[0],gnx[i],index[0],index[i],bGroup,
 373                   &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 374         fprintf(dist,"  %12e",bMin?dmin:dmax);
 375         if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 376         if (nres) {
 377           for(j=0; j<nres; j++) {
 378             calc_dist(rcut,bPBC,ePBC,box,x0,residue[j+1]-residue[j],gnx[i],
 379                       &(index[0][residue[j]]),index[i],bGroup,
 380                       &dmin,&dmax,&nmin,&nmax,&min1r,&min2r,&max1r,&max2r);
 381             mindres[i-1][j] = min(mindres[i-1][j],dmin);
 382             maxdres[i-1][j] = max(maxdres[i-1][j],dmax);
 383           }
 384         }
 385       }
 386     }
 387     fprintf(dist,"\n");
 388     if (num)
 389       fprintf(num,"\n");
 390     if ( bMin?min1:max1 != -1 )
 391       if (atm)
 392         fprintf(atm,"%12e  %12d  %12d\n",
 393                 output_env_conv_time(oenv,t),1+(bMin ? min1 : max1),
 394                                   1+(bMin ? min2 : max2));
 395
 396     if (trxout) {
 397       oindex[0]=bMin?min1:max1;
 398       oindex[1]=bMin?min2:max2;
 399       write_trx(trxout,2,oindex,atoms,i,t,box,x0,NULL,NULL);
 400     }
 401     bFirst=FALSE;
 402         /*dmin should be minimum distance for residue and group*/
 403         if (bEachResEachTime)
 404         {
 405             fprintf(respertime,"%12e",t);
 406             for(i=1; i<ng; i++)
 407                         for(j=0; j<nres; j++)
 408                         {
 409                                 fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 410                                 /*reset distances for next time point*/
 411                                 mindres[i-1][j]=1e6;
 412                                 maxdres[i-1][j]=0;
 413                         }
 414                 fprintf(respertime, "\n");
 415         }
 416   } while (read_next_x(oenv,status,&t,natoms,x0,box));
 417
 418   close_trj(status);
 419   ffclose(dist);
 420   if (num) ffclose(num);
 421   if (atm) ffclose(atm);
 422   if (trxout) close_trx(trxout);
 423
 424   if(nres && !bEachResEachTime) {
 425     FILE *res;
 426
 427     sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
 428     res=xvgropen(rfile,buf,"Residue (#)","Distance (nm)",oenv);
 429     xvgr_legend(res,ng-1,(const char**)leg,oenv);
 430     for(j=0; j<nres; j++) {
 431       fprintf(res, "%4d", j+1);
 432       for(i=1; i<ng; i++) {
 433         fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 434       }
 435       fprintf(res, "\n");
 436     }
 437   }
 438
 439   sfree(x0);
 440 }
 441
 442 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
 443 {
 444   int i;
 445   int nres=0,resnr, presnr;
 446   int *residx;
 447
 448   /* build index of first atom numbers for each residue */
 449   presnr = NOTSET;
 450   snew(residx, atoms->nres+1);
 451   for(i=0; i<n; i++) {
 452     resnr = atoms->atom[index[i]].resind;
 453     if (resnr != presnr) {
 454       residx[nres]=i;
 455       nres++;
 456       presnr=resnr;
 457     }
 458   }
 459   if (debug) printf("Found %d residues out of %d (%d/%d atoms)\n",
 460                     nres, atoms->nres, atoms->nr, n);
 461   srenew(residx, nres+1);
 462   /* mark end of last residue */
 463   residx[nres] = n;
 464   *resindex = residx;
 465   return nres;
 466 }
 467
 468 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
 469 {
 470   int i,j;
 471
 472   for(i=0; i<nres-1; i++) {
 473     fprintf(out,"Res %d (%d):", i, resindex[i+1]-resindex[i]);
 474     for(j=resindex[i]; j<resindex[i+1]; j++)
 475       fprintf(out," %d(%d)", j, index[j]);
 476     fprintf(out,"\n");
 477   }
 478 }
 479
 480 int gmx_mindist(int argc,char *argv[])
 481 {
 482   const char *desc[] = {
 483     "g_mindist computes the distance between one group and a number of",
 484     "other groups. Both the minimum distance",
 485     "(between any pair of atoms from the respective groups)",
 486     "and the number of contacts within a given",
 487     "distance are written to two separate output files.",
 488     "With the [TT]-group[tt] option a contact of an atom an other group",
 489     "with multiple atoms in the first group is counted as one contact",
 490     "instead of as multiple contacts.",
 491     "With [TT]-or[tt], minimum distances to each residue in the first",
 492     "group are determined and plotted as a function of residue number.[PAR]",
 493     "With option [TT]-pi[tt] the minimum distance of a group to its",
 494     "periodic image is plotted. This is useful for checking if a protein",
 495     "has seen its periodic image during a simulation. Only one shift in",
 496     "each direction is considered, giving a total of 26 shifts.",
 497     "It also plots the maximum distance within the group and the lengths",
 498     "of the three box vectors.[PAR]",
 499     "Other programs that calculate distances are [TT]g_dist[tt]",
 500     "and [TT]g_bond[tt]."
 501   };
 502   const char *bugs[] = {
 503     "The [TT]-pi[tt] option is very slow."
 504   };
 505
 506   static gmx_bool bMat=FALSE,bPI=FALSE,bSplit=FALSE,bMax=FALSE,bPBC=TRUE;
 507   static gmx_bool bGroup=FALSE;
 508   static real rcutoff=0.6;
 509   static int  ng=1;
 510   static gmx_bool bEachResEachTime=FALSE,bPrintResName=FALSE;
 511   t_pargs pa[] = {
 512     { "-matrix", FALSE, etBOOL, {&bMat},
 513       "Calculate half a matrix of group-group distances" },
 514     { "-max",    FALSE, etBOOL, {&bMax},
 515       "Calculate *maximum* distance instead of minimum" },
 516     { "-d",      FALSE, etREAL, {&rcutoff},
 517       "Distance for contacts" },
 518     { "-group",      FALSE, etBOOL, {&bGroup},
 519       "Count contacts with multiple atoms in the first group as one" },
 520     { "-pi",     FALSE, etBOOL, {&bPI},
 521       "Calculate minimum distance with periodic images" },
 522     { "-split",  FALSE, etBOOL, {&bSplit},
 523       "Split graph where time is zero" },
 524     { "-ng",       FALSE, etINT, {&ng},
 525       "Number of secondary groups to compute distance to a central group" },
 526     { "-pbc",    FALSE, etBOOL, {&bPBC},
 527       "Take periodic boundary conditions into account" },
 528     { "-respertime",  FALSE, etBOOL, {&bEachResEachTime},
 529       "When writing per-residue distances, write distance for each time point" },
 530     { "-printresname",  FALSE, etBOOL, {&bPrintResName},
 531       "Write residue names" }
 532   };
 533   output_env_t oenv;
 534   t_topology *top=NULL;
 535   int        ePBC=-1;
 536   char       title[256];
 537   real       t;
 538   rvec       *x;
 539   matrix     box;
 540   gmx_bool       bTop=FALSE;
 541
 542   FILE      *atm;
 543   int       i,j,nres=0;
 544   const char *trxfnm,*tpsfnm,*ndxfnm,*distfnm,*numfnm,*atmfnm,*oxfnm,*resfnm;
 545   char      **grpname;
 546   int       *gnx;
 547   atom_id   **index, *residues=NULL;
 548   t_filenm  fnm[] = {
 549     { efTRX, "-f",  NULL,      ffREAD },
 550     { efTPS,  NULL, NULL,      ffOPTRD },
 551     { efNDX,  NULL, NULL,      ffOPTRD },
 552     { efXVG, "-od","mindist",  ffWRITE },
 553     { efXVG, "-on","numcont",  ffOPTWR },
 554     { efOUT, "-o", "atm-pair", ffOPTWR },
 555     { efTRO, "-ox","mindist",  ffOPTWR },
 556     { efXVG, "-or","mindistres", ffOPTWR }
 557   };
 558 #define NFILE asize(fnm)
 559
 560   CopyRight(stderr,argv[0]);
 561   parse_common_args(&argc,argv,
 562                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 563                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL, &oenv);
 564
 565   trxfnm = ftp2fn(efTRX,NFILE,fnm);
 566   ndxfnm = ftp2fn_null(efNDX,NFILE,fnm);
 567   distfnm= opt2fn("-od",NFILE,fnm);
 568   numfnm = opt2fn_null("-on",NFILE,fnm);
 569   atmfnm = ftp2fn_null(efOUT,NFILE,fnm);
 570   oxfnm  = opt2fn_null("-ox",NFILE,fnm);
 571   resfnm = opt2fn_null("-or",NFILE,fnm);
 572   if (bPI || resfnm != NULL) {
 573     /* We need a tps file */
 574     tpsfnm = ftp2fn(efTPS,NFILE,fnm);
 575   } else {
 576     tpsfnm = ftp2fn_null(efTPS,NFILE,fnm);
 577   }
 578
 579   if (!tpsfnm && !ndxfnm)
 580     gmx_fatal(FARGS,"You have to specify either the index file or a tpr file");
 581
 582   if (bPI) {
 583     ng = 1;
 584     fprintf(stderr,"Choose a group for distance calculation\n");
 585   }
 586   else if (!bMat)
 587     ng++;
 588
 589   snew(gnx,ng);
 590   snew(index,ng);
 591   snew(grpname,ng);
 592
 593   if (tpsfnm || resfnm || !ndxfnm) {
 594     snew(top,1);
 595     bTop = read_tps_conf(tpsfnm,title,top,&ePBC,&x,NULL,box,FALSE);
 596     if (bPI && !bTop)
 597       printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
 598   }
 599   get_index(top ? &(top->atoms) : NULL,ndxfnm,ng,gnx,index,grpname);
 600
 601   if (bMat && (ng == 1)) {
 602     ng = gnx[0];
 603     printf("Special case: making distance matrix between all atoms in group %s\n",
 604            grpname[0]);
 605     srenew(gnx,ng);
 606     srenew(index,ng);
 607     srenew(grpname,ng);
 608     for(i=1; (i<ng); i++) {
 609       gnx[i]      = 1;
 610       grpname[i]  = grpname[0];
 611       snew(index[i],1);
 612       index[i][0] = index[0][i];
 613     }
 614     gnx[0] = 1;
 615   }
 616
 617   if (resfnm) {
 618     nres=find_residues(top ? &(top->atoms) : NULL,
 619                        gnx[0], index[0], &residues);
 620     if (debug) dump_res(debug, nres, residues, gnx[0], index[0]);
 621   }
 622
 623   if (bPI) {
 624     periodic_mindist_plot(trxfnm,distfnm,top,ePBC,gnx[0],index[0],bSplit,oenv);
 625   } else {
 626     dist_plot(trxfnm,atmfnm,distfnm,numfnm,resfnm,oxfnm,
 627               rcutoff,bMat,top ? &(top->atoms) : NULL,
 628               ng,index,gnx,grpname,bSplit,!bMax, nres, residues,bPBC,ePBC,
 629               bGroup,bEachResEachTime,bPrintResName,oenv);
 630   }
 631
 632   do_view(oenv,distfnm,"-nxy");
 633   if (!bPI)
 634     do_view(oenv,numfnm,"-nxy");
 635
 636   thanx(stderr);
 637
 638   return 0;
 639 }
 640