src/tools/g_mindist.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33
  34 #include <math.h>
  35 #include <stdlib.h>
  36 #include "sysstuff.h"
  37 #include "string.h"
  38 #include "typedefs.h"
  39 #include "smalloc.h"
  40 #include "macros.h"
  41 #include "vec.h"
  42 #include "xvgr.h"
  43 #include "pbc.h"
  44 #include "copyrite.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "rdgroup.h"
  48 #include "tpxio.h"
  49 #include "rmpbc.h"
  50
  51 static void periodic_dist(matrix box,rvec x[],int n,atom_id index[],
  52                           real *rmin,real *rmax)
  53 {
  54 #define NSHIFT 26
  55   int  sx,sy,sz,i,j,s;
  56   real sqr_box,r2min,r2max,r2;
  57   rvec shift[NSHIFT],d0,d;
  58
  59   sqr_box = sqr(min(box[XX][XX],min(box[YY][YY],box[ZZ][ZZ])));
  60
  61   s = 0;
  62   for(sz=-1; sz<=1; sz++)
  63     for(sy=-1; sy<=1; sy++)
  64       for(sx=-1; sx<=1; sx++)
  65         if (sx!=0 || sy!=0 || sz!=0) {
  66           for(i=0; i<DIM; i++)
  67             shift[s][i] = sx*box[XX][i]+sy*box[YY][i]+sz*box[ZZ][i];
  68           s++;
  69         }
  70
  71   r2min = sqr_box;
  72   r2max = 0;
  73
  74   for(i=0; i<n; i++)
  75     for(j=i+1; j<n; j++) {
  76       rvec_sub(x[index[i]],x[index[j]],d0);
  77       r2 = norm2(d0);
  78       if (r2 > r2max)
  79         r2max = r2;
  80       for(s=0; s<NSHIFT; s++) {
  81         rvec_add(d0,shift[s],d);
  82         r2 = norm2(d);
  83         if (r2 < r2min)
  84           r2min = r2;
  85       }
  86     }
  87
  88   *rmin = sqrt(r2min);
  89   *rmax = sqrt(r2max);
  90 }
  91
  92 static void periodic_mindist_plot(char *trxfn,char *outfn,
  93                                   t_topology *top,int n,atom_id index[],
  94                                   bool bSplit)
  95 {
  96   FILE   *out;
  97   char   *leg[5] = { "min per.","max int.","box1","box2","box3" };
  98   int    status;
  99   real   t;
 100   rvec   *x;
 101   matrix box;
 102   int    natoms;
 103   real   r,rmin,rmax,rmint,tmint;
 104   bool   bFirst;
 105
 106   natoms=read_first_x(&status,trxfn,&t,&x,box);
 107
 108   check_index(NULL,n,index,NULL,natoms);
 109
 110   out = xvgropen(outfn,"Minimum distance to periodic image",
 111                  time_label(),"Distance (nm)");
 112   fprintf(out,"@ subtitle \"and maximum internal distance\"\n");
 113   xvgr_legend(out,5,leg);
 114
 115   rmint = box[XX][XX];
 116   tmint = 0;
 117
 118   bFirst=TRUE;
 119   do {
 120     rm_pbc(&(top->idef),natoms,box,x,x);
 121     periodic_dist(box,x,n,index,&rmin,&rmax);
 122     if (rmin < rmint) {
 123       rmint = rmin;
 124       tmint = t;
 125     }
 126     if ( bSplit && !bFirst && abs(t/time_factor())<1e-5 )
 127       fprintf(out, "&\n");
 128     fprintf(out,"\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
 129             convert_time(t),rmin,rmax,norm(box[0]),norm(box[1]),norm(box[2]));
 130     bFirst=FALSE;
 131   } while(read_next_x(status,&t,natoms,x,box));
 132
 133   fclose(out);
 134
 135   fprintf(stdout,
 136           "\nThe shortest periodic distance is %g (nm) at time %g (%s)\n",
 137           rmint,convert_time(tmint),time_unit());
 138 }
 139
 140 static void calc_dist(real rcut, matrix box, rvec x[],
 141                       int nx1,int nx2, atom_id index1[], atom_id index2[],
 142                       real *rmin, real *rmax, int *nmin, int *nmax,
 143                       int *ixmin, int *jxmin, int *ixmax, int *jxmax)
 144 {
 145   int     i,j,j0=0,j1;
 146   int     ix,jx;
 147   atom_id *index3;
 148   rvec    dx;
 149   real    r2,rmin2,rmax2,rcut2;
 150
 151   *ixmin = -1;
 152   *jxmin = -1;
 153   *ixmax = -1;
 154   *jxmax = -1;
 155   *nmin = 0;
 156   *nmax = 0;
 157
 158   rcut2=sqr(rcut);
 159
 160   /* Must init pbc every step because of pressure coupling */
 161   init_pbc(box);
 162   if (index2) {
 163     j0=0;
 164     j1=nx2;
 165     index3=index2;
 166   } else {
 167     j1=nx1;
 168     index3=index1;
 169   }
 170
 171   rmin2=1e12;
 172   rmax2=-1e12;
 173
 174   for(i=0; (i < nx1); i++) {
 175     ix=index1[i];
 176     if (!index2)
 177       j0=i+1;
 178     for(j=j0; (j < j1); j++) {
 179       jx=index3[j];
 180       if (ix != jx) {
 181         pbc_dx(x[ix],x[jx],dx);
 182         r2=iprod(dx,dx);
 183         if (r2 < rmin2) {
 184           rmin2=r2;
 185           *ixmin=ix;
 186           *jxmin=jx;
 187         }
 188         if (r2 > rmax2) {
 189           rmax2=r2;
 190           *ixmax=ix;
 191           *jxmax=jx;
 192         }
 193         if (r2 < rcut2)
 194           *nmin++;
 195         else if (r2 > rcut2)
 196           *nmax++;
 197       }
 198     }
 199   }
 200   *rmin = sqrt(rmin2);
 201   *rmax = sqrt(rmax2);
 202 }
 203
 204 void dist_plot(char *fn,char *afile,char *dfile,
 205                char *nfile,char *rfile,char *xfile,
 206                real rcut,bool bMat,t_atoms *atoms,
 207                int ng,atom_id *index[],int gnx[],char *grpn[],bool bSplit,
 208                bool bMin, int nres, atom_id *residue)
 209 {
 210   FILE         *atm,*dist,*num;
 211   int          trxout;
 212   char         buf[256];
 213   char         **leg;
 214 <<<<<<< g_mindist.c
 215   real         t,md,**mdist=NULL;
 216   int          nd,status;
 217 =======
 218   real         t,dmin,dmax,**mindres,**maxdres;
 219   int          nmin,nmax,status;
 220 >>>>>>> 1.29
 221   int          i,j,k,natoms;
 222   int          min1,min2,max1,max2;
 223   atom_id      oindex[2];
 224   rvec         *x0;
 225   matrix       box;
 226   t_trxframe   frout;
 227   bool         bFirst;
 228
 229   if ((natoms=read_first_x(&status,fn,&t,&x0,box))==0)
 230     fatal_error(0,"Could not read coordinates from statusfile\n");
 231
 232   sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
 233   dist= xvgropen(dfile,buf,time_label(),"Distance (nm)");
 234   sprintf(buf,"Number of Contacts %s %g nm",bMin ? "<" : ">",rcut);
 235   num = nfile ? xvgropen(nfile,buf,time_label(),"Number") : NULL;
 236   atm = afile ? ffopen(afile,"w") : NULL;
 237   trxout = xfile ? open_trx(xfile,"w") : NOTSET;
 238
 239   if (bMat) {
 240     if (ng == 1) {
 241       snew(leg,1);
 242       sprintf(buf,"Internal in %s",grpn[0]);
 243       leg[0]=strdup(buf);
 244       xvgr_legend(dist,0,leg);
 245       if (num) xvgr_legend(num,0,leg);
 246     }
 247     else {
 248       snew(leg,(ng*(ng-1))/2);
 249       for(i=j=0; (i<ng-1); i++) {
 250         for(k=i+1; (k<ng); k++,j++) {
 251           sprintf(buf,"%s-%s",grpn[i],grpn[k]);
 252           leg[j]=strdup(buf);
 253         }
 254       }
 255       xvgr_legend(dist,j,leg);
 256       if (num) xvgr_legend(num,j,leg);
 257     }
 258   }
 259   else {
 260     snew(leg,ng-1);
 261     for(i=0; (i<ng-1); i++) {
 262       sprintf(buf,"%s-%s",grpn[0],grpn[i+1]);
 263       leg[i]=strdup(buf);
 264     }
 265     xvgr_legend(dist,ng-1,leg);
 266     if (num) xvgr_legend(num,ng-1,leg);
 267   }
 268   j=0;
 269   if (nres) {
 270     snew(mindres, ng-1);
 271     snew(maxdres, ng-1);
 272     for(i=1; i<ng; i++) {
 273       snew(mindres[i-1], nres);
 274       snew(maxdres[i-1], nres);
 275       for(j=0; j<nres; j++)
 276         mindres[i-1][j]=1e6;
 277       /* maxdres[*][*] is already 0 */
 278     }
 279   }
 280   bFirst=TRUE;
 281   do {
 282     if ( bSplit && !bFirst && abs(t/time_factor())<1e-5 ) {
 283       fprintf(dist, "&\n");
 284       if (num) fprintf(num, "&\n");
 285       if (atm) fprintf(atm, "&\n");
 286     }
 287     fprintf(dist,"%12g",convert_time(t));
 288     if (num) fprintf(num,"%12g",convert_time(t));
 289
 290     if (bMat) {
 291       if (ng == 1) {
 292         calc_dist(rcut,box,x0,gnx[0],gnx[0],index[0],index[0],
 293                   &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 294         fprintf(dist,"  %12g",bMin?dmin:dmax);
 295         if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 296       }
 297       else {
 298         for(i=0; (i<ng-1); i++) {
 299           for(k=i+1; (k<ng); k++) {
 300             calc_dist(rcut,box,x0,gnx[i],gnx[k],index[i],index[k],
 301                       &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 302             fprintf(dist,"  %12g",bMin?dmin:dmax);
 303             if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 304           }
 305         }
 306       }
 307     }
 308     else {
 309       for(i=1; (i<ng); i++) {
 310         calc_dist(rcut,box,x0,gnx[0],gnx[i],index[0],index[i],
 311                   &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 312         fprintf(dist,"  %12g",bMin?dmin:dmax);
 313         if (num) fprintf(num,"  %8d",bMin?nmin:nmax);
 314         if (nres) {
 315           for(j=0; j<nres; j++) {
 316             calc_dist(rcut,box,x0,residue[j+1]-residue[j],gnx[i],
 317                       &(index[0][residue[j]]),index[i],
 318                       &dmin,&dmax,&nmin,&nmax,&min1,&min2,&max1,&max2);
 319             mindres[i-1][j] = min(mindres[i-1][j],dmin);
 320             maxdres[i-1][j] = max(maxdres[i-1][j],dmax);
 321           }
 322         }
 323       }
 324     }
 325     fprintf(dist,"\n");
 326     if (num)
 327       fprintf(num,"\n");
 328     if ( bMin?min1:max1 != -1 )
 329       if (atm)
 330         fprintf(atm,"%12g  %12d  %12d\n",
 331                 convert_time(t),bMin?min1:max1+1,bMin?min2:max2+1);
 332
 333     if (trxout>=0) {
 334       oindex[0]=bMin?min1:max1;
 335       oindex[1]=bMin?min2:max2;
 336       write_trx(trxout,2,oindex,atoms,i,t,box,x0,NULL);
 337     }
 338     bFirst=FALSE;
 339   } while (read_next_x(status,&t,natoms,x0,box));
 340
 341   close_trj(status);
 342   fclose(dist);
 343   if (num) fclose(num);
 344   if (atm) fclose(atm);
 345   if (trxout>=0) close_xtc(trxout);
 346
 347   if(nres) {
 348     FILE *res;
 349
 350     sprintf(buf,"%simum Distance",bMin ? "Min" : "Max");
 351     res=xvgropen(rfile,buf,"Residue (#)","Distance (nm)");
 352     xvgr_legend(res,ng-1,leg);
 353     for(j=0; j<nres; j++) {
 354       fprintf(res, "%4d", j+1);
 355       for(i=1; i<ng; i++) {
 356         fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
 357       }
 358       fprintf(res, "\n");
 359     }
 360   }
 361
 362   sfree(x0);
 363 }
 364
 365 int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
 366 {
 367   int i;
 368   int nres=0,resnr, presnr;
 369   int *residx;
 370
 371   /* build index of first atom numbers for each residue */
 372   presnr = NOTSET;
 373   snew(residx, atoms->nres);
 374   for(i=0; i<n; i++) {
 375     resnr = atoms->atom[index[i]].resnr;
 376     if (resnr != presnr) {
 377       residx[nres]=i;
 378       nres++;
 379       presnr=resnr;
 380     }
 381   }
 382   if (debug) printf("Found %d residues out of %d (%d/%d atoms)\n",
 383                     nres, atoms->nres, atoms->nr, n);
 384   srenew(residx, nres+1);
 385   /* mark end of last residue */
 386   residx[nres]=n+1;
 387   *resindex = residx;
 388   return nres;
 389 }
 390
 391 void dump_res(FILE *out, int nres, atom_id *resindex, int n, atom_id index[])
 392 {
 393   int i,j;
 394
 395   for(i=0; i<nres-1; i++) {
 396     fprintf(out,"Res %d (%d):", i, resindex[i+1]-resindex[i]);
 397     for(j=resindex[i]; j<resindex[i+1]; j++)
 398       fprintf(out," %d(%d)", j, index[j]);
 399     fprintf(out,"\n");
 400   }
 401 }
 402
 403 int main(int argc,char *argv[])
 404 {
 405   static char *desc[] = {
 406     "g_mindist computes the distance between one group and a number of",
 407     "other groups.",
 408     "Both the minimum distance and the number of contacts within a given",
 409     "distance are written to two separate output files.",
 410     "With [TT]-or[tt], minimum distances to each residue in the first",
 411     "group are determined and plotted as a function of reisdue number.[PAR]",
 412     "With option [TT]-pi[tt] the minimum distance of a group to its",
 413     "periodic image is plotted. This is useful for checking if a protein",
 414     "has seen its periodic image during a simulation. Only one shift in",
 415     "each direction is considered, giving a total of 26 shifts.",
 416     "It also plots the maximum distance within the group and the lengths",
 417     "of the three box vectors.[PAR]",
 418     "Other programs that calculate distances are [TT]g_dist[tt]",
 419     "and [TT]g_bond[tt]."
 420   };
 421   static char *bugs[] = {
 422     "The [TT]-pi[tt] option is very slow."
 423   };
 424
 425   static bool bMat=FALSE,bPer=FALSE,bSplit=FALSE,bMax=FALSE;
 426   static real rcutoff=0.6;
 427   t_pargs pa[] = {
 428     { "-matrix", FALSE, etBOOL, {&bMat},
 429       "Calculate half a matrix of group-group distances" },
 430     { "-max",    FALSE, etBOOL, {&bMax},
 431       "Calculate *maximum* distance instead of minimum" },
 432     { "-d",      FALSE, etREAL, {&rcutoff},
 433       "Distance for contacts" },
 434     { "-pi",     FALSE, etBOOL, {&bPer},
 435       "Calculate minimum distance with periodic images" },
 436     { "-split",  FALSE, etBOOL, {&bSplit},
 437       "Split graph where time is zero" },
 438   };
 439   t_topology top;
 440   char       title[256];
 441   real       t;
 442   rvec       *x;
 443   matrix     box;
 444
 445   FILE      *atm;
 446   int       i,j,ng,nres=0;
 447   char      *trxfnm,*tpsfnm,*ndxfnm,*distfnm,*numfnm,*atmfnm,*oxfnm,*resfnm;
 448   char      **grpname;
 449   int       *gnx;
 450   atom_id   **index, *residues=NULL;
 451   t_filenm  fnm[] = {
 452     { efTRX, "-f",  NULL,      ffREAD },
 453     { efTPS,  NULL, NULL,      ffOPTRD },
 454     { efNDX,  NULL, NULL,      ffOPTRD },
 455     { efXVG, "-od","mindist",  ffWRITE },
 456     { efXVG, "-on","numcont",  ffOPTWR },
 457     { efOUT, "-o", "atm-pair", ffOPTWR },
 458     { efTRX, "-ox","mindist",  ffOPTWR },
 459     { efXVG, "-or","mindistres", ffOPTWR }
 460   };
 461 #define NFILE asize(fnm)
 462
 463   CopyRight(stderr,argv[0]);
 464   parse_common_args(&argc,argv,
 465                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
 466                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 467
 468   trxfnm = ftp2fn(efTRX,NFILE,fnm);
 469   tpsfnm = ftp2fn_null(efTPS,NFILE,fnm);
 470   ndxfnm = ftp2fn_null(efNDX,NFILE,fnm);
 471   distfnm= opt2fn("-od",NFILE,fnm);
 472   numfnm = opt2fn_null("-on",NFILE,fnm);
 473   atmfnm = ftp2fn_null(efOUT,NFILE,fnm);
 474   oxfnm  = opt2fn_null("-ox",NFILE,fnm);
 475   resfnm = opt2fn_null("-or",NFILE,fnm);
 476
 477   if (bPer) {
 478     ng = 1;
 479     fprintf(stderr,"Choose a group for distance calculation\n");
 480   } else {
 481     if (bMat)
 482       fprintf(stderr,"You can compute all distances between a number of groups\n"
 483               "How many groups do you want (>= 1) ?\n");
 484     else
 485       fprintf(stderr,"You can compute the distances between a first group\n"
 486               "and a number of other groups.\n"
 487               "How many other groups do you want (>= 1) ?\n");
 488     ng = 0;
 489     do {
 490       scanf("%d",&ng);
 491       if (!bMat)
 492         ng++;
 493     } while (ng < 1);
 494   }
 495   snew(gnx,ng);
 496   snew(index,ng);
 497   snew(grpname,ng);
 498
 499   if (tpsfnm || !ndxfnm)
 500     read_tps_conf(tpsfnm,title,&top,&x,NULL,box,FALSE);
 501
 502   get_index(&top.atoms,ndxfnm,ng,gnx,index,grpname);
 503
 504   if (bMat && (ng == 1)) {
 505     ng = gnx[0];
 506     printf("Special case: making distance matrix between all atoms in group %s\n",
 507            grpname[0]);
 508     srenew(gnx,ng);
 509     srenew(index,ng);
 510     srenew(grpname,ng);
 511     for(i=1; (i<ng); i++) {
 512       gnx[i]      = 1;
 513       grpname[i]  = grpname[0];
 514       snew(index[i],1);
 515       index[i][0] = index[0][i];
 516     }
 517     gnx[0] = 1;
 518   }
 519
 520   if (resfnm) {
 521     nres=find_residues(&top.atoms, gnx[0], index[0], &residues);
 522     if (debug) dump_res(debug, nres, residues, gnx[0], index[0]);
 523   }
 524
 525   if (bPer)
 526     periodic_mindist_plot(trxfnm,distfnm,&top,gnx[0],index[0],bSplit);
 527   else
 528     dist_plot(trxfnm,atmfnm,distfnm,numfnm,resfnm,oxfnm,
 529               rcutoff,bMat,&top.atoms,ng,index,gnx,grpname,bSplit,
 530               !bMax, nres, residues);
 531
 532   do_view(distfnm,"-nxy");
 533   if (!bPer)
 534     do_view(numfnm,"-nxy");
 535
 536   thanx(stderr);
 537
 538   return 0;
 539 }
 540