src/tools/g_mdmat.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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_mdmat_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "macros.h"
  34 #include "vec.h"
  35 #include "sysstuff.h"
  36 #include "typedefs.h"
  37 #include "filenm.h"
  38 #include "statutil.h"
  39 #include "copyrite.h"
  40 #include "futil.h"
  41 #include "fatal.h"
  42 #include "smalloc.h"
  43 #include "matio.h"
  44 #include "xvgr.h"
  45 #include "rdgroup.h"
  46 #include "tpxio.h"
  47 #include "rmpbc.h"
  48
  49 #define FARAWAY 10000
  50
  51 int *res_ndx(t_atoms *atoms)
  52 {
  53   int *rndx;
  54   int i,r0;
  55
  56   if (atoms->nr <= 0)
  57     return NULL;
  58   snew(rndx,atoms->nr);
  59   r0=atoms->atom[0].resnr;
  60   for(i=0; (i<atoms->nr); i++)
  61     rndx[i]=atoms->atom[i].resnr-r0;
  62
  63   return rndx;
  64 }
  65
  66 int *res_natm(t_atoms *atoms)
  67 {
  68   int *natm;
  69   int i,j,r0;
  70
  71   if (atoms->nr <= 0)
  72     return NULL;
  73   snew(natm,atoms->nres);
  74   r0=atoms->atom[0].resnr;
  75   j=0;
  76   for(i=0; (i<atoms->nres); i++) {
  77     while ((atoms->atom[j].resnr)-r0 == i) {
  78       natm[i]++;
  79       j++;
  80     }
  81   }
  82
  83   return natm;
  84 }
  85
  86 static void calc_mat(int nres, int natoms, int rndx[],
  87                      rvec x[], atom_id *index,
  88                      real trunc, real **mdmat,int **nmat)
  89 {
  90   int i,j,resi,resj;
  91   real trunc2,r,r2;
  92
  93   trunc2=sqr(trunc);
  94   for(resi=0; (resi<nres); resi++)
  95     for(resj=0; (resj<nres); resj++)
  96       mdmat[resi][resj]=FARAWAY;
  97   for(i=0; (i<natoms); i++) {
  98     resi=rndx[i];
  99     for(j=i+1; (j<natoms); j++) {
 100       resj=rndx[j];
 101       r2=distance2(x[index[i]],x[index[j]]);
 102       if ( r2 < trunc2 ) {
 103         nmat[resi][j]++;
 104         nmat[resj][i]++;
 105       }
 106       mdmat[resi][resj]=min(r2,mdmat[resi][resj]);
 107     }
 108   }
 109
 110   for(resi=0; (resi<nres); resi++) {
 111     mdmat[resi][resi]=0;
 112     for(resj=resi+1; (resj<nres); resj++) {
 113       r=sqrt(mdmat[resi][resj]);
 114       mdmat[resi][resj]=r;
 115       mdmat[resj][resi]=r;
 116     }
 117   }
 118 }
 119
 120 static void tot_nmat(int nres, int natoms, int nframes, int **nmat,
 121                      int *tot_n, real *mean_n)
 122 {
 123   int i,j;
 124
 125   for (i=0; (i<nres); i++) {
 126     for (j=0; (j<natoms); j++)
 127       if (nmat[i][j] != 0) {
 128         tot_n[i]++;
 129         mean_n[i]+=nmat[i][j];
 130       }
 131     mean_n[i]/=nframes;
 132   }
 133 }
 134
 135 int main(int argc,char *argv[])
 136 {
 137   static char *desc[] = {
 138     "g_mdmat makes distance matrices consisting of the smallest distance",
 139     "between residue pairs. With -frames these distance matrices can be",
 140     "stored as a function",
 141     "of time, to be able to see differences in tertiary structure as a",
 142     "funcion of time. If you choose your options unwise, this may generate",
 143     "a large output file. Default only an averaged matrix over the whole",
 144     "trajectory is output.",
 145     "Also a count of the number of different atomic contacts between",
 146     "residues over the whole trajectory can be made.",
 147     "The output can be processed with xpm2ps to make a PostScript (tm) plot."
 148   };
 149   static real truncate=1.5,dt=0.0;
 150   static bool bAtom=FALSE;
 151   static int  nlevels=40;
 152   t_pargs pa[] = {
 153     { "-t",   FALSE, etREAL, {&truncate},
 154       "trunc distance" },
 155     { "-nlevels",   FALSE, etINT,  {&nlevels},
 156       "Discretize distance in # levels" },
 157     { "-dt",  FALSE, etREAL, {&dt},
 158       "Only analyze a frame each dt picoseconds" }
 159   };
 160   t_filenm   fnm[] = {
 161     { efTRX, "-f",  NULL, ffREAD },
 162     { efTPS, NULL,  NULL, ffREAD },
 163     { efNDX, NULL,  NULL, ffOPTRD },
 164     { efXPM, "-mean", "dm", ffWRITE },
 165     { efXPM, "-frames", "dmf", ffOPTWR },
 166     { efXVG, "-no", "num",ffOPTWR },
 167   };
 168 #define NFILE asize(fnm)
 169
 170   FILE       *out=NULL,*fp;
 171   t_topology top;
 172   t_atoms    useatoms;
 173   int        isize;
 174   atom_id    *index;
 175   char       *grpname;
 176   int        *rndx,*natm;
 177
 178   int        i,j,status,nres,natoms,nframes,it,teller,trxnat;
 179   int        nr0;
 180   bool       bCalcN,bFrames;
 181   real       t,t0,ratio;
 182   char       title[256],label[234];
 183   t_rgb      rlo,rhi;
 184   rvec       *x;
 185   real       **mdmat,*resnr,**totmdmat;
 186   int        **nmat,**totnmat;
 187   real       *mean_n;
 188   int        *tot_n;
 189   matrix     box;
 190
 191   CopyRight(stdout,argv[0]);
 192
 193   parse_common_args(&argc,argv,PCA_CAN_TIME,TRUE,NFILE,fnm,
 194                     asize(pa),pa,asize(desc),desc,0,NULL);
 195
 196   fprintf(stderr,"Will truncate at %f nm\n",truncate);
 197   bCalcN = opt2bSet("-no",NFILE,fnm);
 198   bFrames= opt2bSet("-frames",NFILE,fnm);
 199   if ( bCalcN )
 200     fprintf(stderr,"Will calculate number of different contacts\n");
 201   fprintf(stderr,"Time interval between frames is %g ps\n",dt);
 202
 203   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&x,NULL,box,FALSE);
 204
 205   fprintf(stderr,"Select group for analysis\n");
 206   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 207
 208   natoms=isize;
 209   snew(useatoms.atom,natoms);
 210   snew(useatoms.atomname,natoms);
 211
 212   useatoms.nres=0;
 213   useatoms.resname=top.atoms.resname;
 214   for(i=0;(i<isize);i++) {
 215     useatoms.atomname[i]=top.atoms.atomname[index[i]];
 216     useatoms.atom[i].resnr=top.atoms.atom[index[i]].resnr;
 217     useatoms.nres=max(useatoms.nres,useatoms.atom[i].resnr+1);
 218   }
 219   useatoms.nr=isize;
 220
 221   rndx=res_ndx(&(useatoms));
 222   natm=res_natm(&(useatoms));
 223   nres=useatoms.nres;
 224   fprintf(stderr,"There are %d residues with %d atoms\n",nres,natoms);
 225
 226   snew(resnr,nres);
 227   snew(mdmat,nres);
 228   snew(nmat,nres);
 229   snew(totnmat,nres);
 230   snew(mean_n,nres);
 231   snew(tot_n,nres);
 232   for(i=0; (i<nres); i++) {
 233     snew(mdmat[i],nres);
 234     snew(nmat[i],natoms);
 235     snew(totnmat[i],natoms);
 236     resnr[i]=i+1;
 237   }
 238   snew(totmdmat,nres);
 239   for(i=0; (i<nres); i++)
 240     snew(totmdmat[i],nres);
 241
 242   trxnat=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 243
 244   teller=0;
 245   nframes=0;
 246
 247   rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
 248   rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
 249   if (bFrames)
 250     out=opt2FILE("-frames",NFILE,fnm,"w");
 251   t0=t;
 252   do {
 253     teller++;
 254     rm_pbc(&top.idef,trxnat,box,x,x);
 255     if (t >= t0) {
 256       nframes++;
 257       calc_mat(nres,natoms,rndx,x,index,truncate,mdmat,nmat);
 258       for (i=0; (i<nres); i++)
 259         for (j=0; (j<natoms); j++)
 260           if (nmat[i][j])
 261             totnmat[i][j]++;
 262       for (i=0; (i<nres); i++)
 263         for (j=0; (j<nres); j++)
 264           totmdmat[i][j] += mdmat[i][j];
 265       if (bFrames) {
 266         sprintf(label,"t=%.0f ps",t);
 267         write_xpm(out,label,"Distance (nm)","Residue Index","Residue Index",
 268                   nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 269       }
 270       t0+=dt;
 271     }
 272   } while (read_next_x(status,&t,trxnat,x,box));
 273   fprintf(stderr,"\n");
 274   close_trj(status);
 275   if (bFrames)
 276     fclose(out);
 277
 278   fprintf(stderr,"Processed %d frames out of %d\n",nframes,teller);
 279
 280   for (i=0; (i<nres); i++)
 281     for (j=0; (j<nres); j++)
 282       totmdmat[i][j] /= nframes;
 283   write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),"Mean smallest distance",
 284             "Distance (nm)","Residue Index","Residue Index",
 285             nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 286
 287   if ( bCalcN ) {
 288     tot_nmat(nres,natoms,nframes,totnmat,tot_n,mean_n);
 289     fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 290                 "Increase in number of contacts","Residue","Ratio");
 291     fprintf(fp,"@ legend on\n");
 292     fprintf(fp,"@ legend box on\n");
 293     fprintf(fp,"@ legend loctype view\n");
 294     fprintf(fp,"@ legend 0.75, 0.8\n");
 295     fprintf(fp,"@ legend string 0 \"Total/mean\"\n");
 296     fprintf(fp,"@ legend string 1 \"Total\"\n");
 297     fprintf(fp,"@ legend string 2 \"Mean\"\n");
 298     fprintf(fp,"@ legend string 3 \"# atoms\"\n");
 299     fprintf(fp,"@ legend string 4 \"Mean/# atoms\"\n");
 300     fprintf(fp,"#%3s %8s  %3s  %8s  %3s  %8s\n",
 301             "res","ratio","tot","mean","natm","mean/atm");
 302     for (i=0; (i<nres); i++) {
 303       if (mean_n[i]==0)
 304         ratio=1;
 305       else
 306         ratio=tot_n[i]/mean_n[i];
 307       fprintf(fp,"%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 308               i+1,ratio,tot_n[i],mean_n[i],natm[i],mean_n[i]/natm[i]);
 309     }
 310     fclose(fp);
 311   }
 312
 313   thanx(stdout);
 314
 315   return 0;
 316 }