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