src/tools/gmx_mdmat.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 <string.h>
  41
  42 #include "macros.h"
  43 #include "vec.h"
  44 #include "sysstuff.h"
  45 #include "typedefs.h"
  46 #include "filenm.h"
  47 #include "statutil.h"
  48 #include "copyrite.h"
  49 #include "futil.h"
  50 #include "gmx_fatal.h"
  51 #include "smalloc.h"
  52 #include "matio.h"
  53 #include "xvgr.h"
  54 #include "index.h"
  55 #include "tpxio.h"
  56 #include "rmpbc.h"
  57 #include "pbc.h"
  58 #include "gmx_ana.h"
  59
  60
  61 #define FARAWAY 10000
  62
  63 int *res_ndx(t_atoms *atoms)
  64 {
  65   int *rndx;
  66   int i,r0;
  67
  68   if (atoms->nr <= 0)
  69     return NULL;
  70   snew(rndx,atoms->nr);
  71   r0=atoms->atom[0].resind;
  72   for(i=0; (i<atoms->nr); i++)
  73     rndx[i]=atoms->atom[i].resind-r0;
  74
  75   return rndx;
  76 }
  77
  78 int *res_natm(t_atoms *atoms)
  79 {
  80   int *natm;
  81   int i,j,r0;
  82
  83   if (atoms->nr <= 0)
  84     return NULL;
  85   snew(natm,atoms->nres);
  86   r0=atoms->atom[0].resind;
  87   j=0;
  88   for(i=0; (i<atoms->nres); i++) {
  89     while ((atoms->atom[j].resind)-r0 == i) {
  90       natm[i]++;
  91       j++;
  92     }
  93   }
  94
  95   return natm;
  96 }
  97
  98 static void calc_mat(int nres, int natoms, int rndx[],
  99                      rvec x[], atom_id *index,
 100                      real trunc, real **mdmat,int **nmat,int ePBC,matrix box)
 101 {
 102   int i,j,resi,resj;
 103   real trunc2,r,r2;
 104   t_pbc pbc;
 105   rvec ddx;
 106
 107   set_pbc(&pbc,ePBC,box);
 108   trunc2=sqr(trunc);
 109   for(resi=0; (resi<nres); resi++)
 110     for(resj=0; (resj<nres); resj++)
 111       mdmat[resi][resj]=FARAWAY;
 112   for(i=0; (i<natoms); i++) {
 113     resi=rndx[i];
 114     for(j=i+1; (j<natoms); j++) {
 115       resj=rndx[j];
 116       pbc_dx(&pbc,x[index[i]],x[index[j]],ddx);
 117       r2 = norm2(ddx);
 118       if ( r2 < trunc2 ) {
 119         nmat[resi][j]++;
 120         nmat[resj][i]++;
 121       }
 122       mdmat[resi][resj]=min(r2,mdmat[resi][resj]);
 123     }
 124   }
 125
 126   for(resi=0; (resi<nres); resi++) {
 127     mdmat[resi][resi]=0;
 128     for(resj=resi+1; (resj<nres); resj++) {
 129       r=sqrt(mdmat[resi][resj]);
 130       mdmat[resi][resj]=r;
 131       mdmat[resj][resi]=r;
 132     }
 133   }
 134 }
 135
 136 static void tot_nmat(int nres, int natoms, int nframes, int **nmat,
 137                      int *tot_n, real *mean_n)
 138 {
 139   int i,j;
 140
 141   for (i=0; (i<nres); i++) {
 142     for (j=0; (j<natoms); j++)
 143       if (nmat[i][j] != 0) {
 144         tot_n[i]++;
 145         mean_n[i]+=nmat[i][j];
 146       }
 147     mean_n[i]/=nframes;
 148   }
 149 }
 150
 151 int gmx_mdmat(int argc,char *argv[])
 152 {
 153   const char *desc[] = {
 154     "g_mdmat makes distance matrices consisting of the smallest distance",
 155     "between residue pairs. With -frames these distance matrices can be",
 156     "stored as a function",
 157     "of time, to be able to see differences in tertiary structure as a",
 158     "funcion of time. If you choose your options unwise, this may generate",
 159     "a large output file. Default only an averaged matrix over the whole",
 160     "trajectory is output.",
 161     "Also a count of the number of different atomic contacts between",
 162     "residues over the whole trajectory can be made.",
 163     "The output can be processed with xpm2ps to make a PostScript (tm) plot."
 164   };
 165   static real truncate=1.5;
 166   static bool bAtom=FALSE;
 167   static int  nlevels=40;
 168   t_pargs pa[] = {
 169     { "-t",   FALSE, etREAL, {&truncate},
 170       "trunc distance" },
 171     { "-nlevels",   FALSE, etINT,  {&nlevels},
 172       "Discretize distance in # levels" }
 173   };
 174   t_filenm   fnm[] = {
 175     { efTRX, "-f",  NULL, ffREAD },
 176     { efTPS, NULL,  NULL, ffREAD },
 177     { efNDX, NULL,  NULL, ffOPTRD },
 178     { efXPM, "-mean", "dm", ffWRITE },
 179     { efXPM, "-frames", "dmf", ffOPTWR },
 180     { efXVG, "-no", "num",ffOPTWR },
 181   };
 182 #define NFILE asize(fnm)
 183
 184   FILE       *out=NULL,*fp;
 185   t_topology top;
 186   int        ePBC;
 187   t_atoms    useatoms;
 188   int        isize;
 189   atom_id    *index;
 190   char       *grpname;
 191   int        *rndx,*natm,prevres,newres;
 192
 193   int        i,j,nres,natoms,nframes,it,trxnat;
 194   t_trxstatus *status;
 195   int        nr0;
 196   bool       bCalcN,bFrames;
 197   real       t,ratio;
 198   char       title[256],label[234];
 199   t_rgb      rlo,rhi;
 200   rvec       *x;
 201   real       **mdmat,*resnr,**totmdmat;
 202   int        **nmat,**totnmat;
 203   real       *mean_n;
 204   int        *tot_n;
 205   matrix     box;
 206   output_env_t oenv;
 207   gmx_rmpbc_t  gpbc=NULL;
 208
 209   CopyRight(stderr,argv[0]);
 210
 211   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,
 212                     asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 213
 214   fprintf(stderr,"Will truncate at %f nm\n",truncate);
 215   bCalcN = opt2bSet("-no",NFILE,fnm);
 216   bFrames= opt2bSet("-frames",NFILE,fnm);
 217   if ( bCalcN )
 218     fprintf(stderr,"Will calculate number of different contacts\n");
 219
 220   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,FALSE);
 221
 222   fprintf(stderr,"Select group for analysis\n");
 223   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 224
 225   natoms=isize;
 226   snew(useatoms.atom,natoms);
 227   snew(useatoms.atomname,natoms);
 228
 229   useatoms.nres = 0;
 230   snew(useatoms.resinfo,natoms);
 231
 232   prevres = top.atoms.atom[index[0]].resind;
 233   newres  = 0;
 234   for(i=0;(i<isize);i++) {
 235     int ii = index[i];
 236     useatoms.atomname[i]=top.atoms.atomname[ii];
 237     if (top.atoms.atom[ii].resind != prevres) {
 238       prevres = top.atoms.atom[ii].resind;
 239       newres++;
 240       useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 241       if (debug) {
 242         fprintf(debug,"New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 243                 *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 244                 *(top.atoms.atomname[ii]),
 245                 ii,i,newres);
 246       }
 247     }
 248     useatoms.atom[i].resind = newres;
 249   }
 250   useatoms.nres = newres+1;
 251   useatoms.nr = isize;
 252
 253   rndx=res_ndx(&(useatoms));
 254   natm=res_natm(&(useatoms));
 255   nres=useatoms.nres;
 256   fprintf(stderr,"There are %d residues with %d atoms\n",nres,natoms);
 257
 258   snew(resnr,nres);
 259   snew(mdmat,nres);
 260   snew(nmat,nres);
 261   snew(totnmat,nres);
 262   snew(mean_n,nres);
 263   snew(tot_n,nres);
 264   for(i=0; (i<nres); i++) {
 265     snew(mdmat[i],nres);
 266     snew(nmat[i],natoms);
 267     snew(totnmat[i],natoms);
 268     resnr[i]=i+1;
 269   }
 270   snew(totmdmat,nres);
 271   for(i=0; (i<nres); i++)
 272     snew(totmdmat[i],nres);
 273
 274   trxnat=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 275
 276   nframes=0;
 277
 278   rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
 279   rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
 280
 281   gpbc = gmx_rmpbc_init(&top.idef,ePBC,trxnat,box);
 282
 283   if (bFrames)
 284     out=opt2FILE("-frames",NFILE,fnm,"w");
 285   do {
 286     gmx_rmpbc(gpbc,trxnat,box,x);
 287     nframes++;
 288     calc_mat(nres,natoms,rndx,x,index,truncate,mdmat,nmat,ePBC,box);
 289     for (i=0; (i<nres); i++)
 290       for (j=0; (j<natoms); j++)
 291         if (nmat[i][j])
 292           totnmat[i][j]++;
 293     for (i=0; (i<nres); i++)
 294       for (j=0; (j<nres); j++)
 295         totmdmat[i][j] += mdmat[i][j];
 296     if (bFrames) {
 297       sprintf(label,"t=%.0f ps",t);
 298       write_xpm(out,0,label,"Distance (nm)","Residue Index","Residue Index",
 299                 nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 300     }
 301   } while (read_next_x(oenv,status,&t,trxnat,x,box));
 302   fprintf(stderr,"\n");
 303   close_trj(status);
 304   gmx_rmpbc_done(gpbc);
 305   if (bFrames)
 306     ffclose(out);
 307
 308   fprintf(stderr,"Processed %d frames\n",nframes);
 309
 310   for (i=0; (i<nres); i++)
 311     for (j=0; (j<nres); j++)
 312       totmdmat[i][j] /= nframes;
 313   write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),0,"Mean smallest distance",
 314             "Distance (nm)","Residue Index","Residue Index",
 315             nres,nres,resnr,resnr,totmdmat,0,truncate,rlo,rhi,&nlevels);
 316
 317   if ( bCalcN ) {
 318     tot_nmat(nres,natoms,nframes,totnmat,tot_n,mean_n);
 319     fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 320                 "Increase in number of contacts","Residue","Ratio",oenv);
 321     fprintf(fp,"@ legend on\n");
 322     fprintf(fp,"@ legend box on\n");
 323     fprintf(fp,"@ legend loctype view\n");
 324     fprintf(fp,"@ legend 0.75, 0.8\n");
 325     fprintf(fp,"@ legend string 0 \"Total/mean\"\n");
 326     fprintf(fp,"@ legend string 1 \"Total\"\n");
 327     fprintf(fp,"@ legend string 2 \"Mean\"\n");
 328     fprintf(fp,"@ legend string 3 \"# atoms\"\n");
 329     fprintf(fp,"@ legend string 4 \"Mean/# atoms\"\n");
 330     fprintf(fp,"#%3s %8s  %3s  %8s  %3s  %8s\n",
 331             "res","ratio","tot","mean","natm","mean/atm");
 332     for (i=0; (i<nres); i++) {
 333       if (mean_n[i]==0)
 334         ratio=1;
 335       else
 336         ratio=tot_n[i]/mean_n[i];
 337       fprintf(fp,"%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 338               i+1,ratio,tot_n[i],mean_n[i],natm[i],mean_n[i]/natm[i]);
 339     }
 340     ffclose(fp);
 341   }
 342
 343   thanx(stderr);
 344
 345   return 0;
 346 }