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 [TT]-frames[tt], these distance matrices can be",
 156     "stored in order to see differences in tertiary structure as a",
 157     "function of time. If you choose your options unwisely, this may generate",
 158     "a large output file. By default, only an averaged matrix over the whole",
 159     "trajectory is output.",
 160     "Also a count of the number of different atomic contacts between",
 161     "residues over the whole trajectory can be made.",
 162     "The output can be processed with xpm2ps to make a PostScript (tm) plot."
 163   };
 164   static real truncate=1.5;
 165   static gmx_bool bAtom=FALSE;
 166   static int  nlevels=40;
 167   t_pargs pa[] = {
 168     { "-t",   FALSE, etREAL, {&truncate},
 169       "trunc distance" },
 170     { "-nlevels",   FALSE, etINT,  {&nlevels},
 171       "Discretize distance in # levels" }
 172   };
 173   t_filenm   fnm[] = {
 174     { efTRX, "-f",  NULL, ffREAD },
 175     { efTPS, NULL,  NULL, ffREAD },
 176     { efNDX, NULL,  NULL, ffOPTRD },
 177     { efXPM, "-mean", "dm", ffWRITE },
 178     { efXPM, "-frames", "dmf", ffOPTWR },
 179     { efXVG, "-no", "num",ffOPTWR },
 180   };
 181 #define NFILE asize(fnm)
 182
 183   FILE       *out=NULL,*fp;
 184   t_topology top;
 185   int        ePBC;
 186   t_atoms    useatoms;
 187   int        isize;
 188   atom_id    *index;
 189   char       *grpname;
 190   int        *rndx,*natm,prevres,newres;
 191
 192   int        i,j,nres,natoms,nframes,it,trxnat;
 193   t_trxstatus *status;
 194   int        nr0;
 195   gmx_bool       bCalcN,bFrames;
 196   real       t,ratio;
 197   char       title[256],label[234];
 198   t_rgb      rlo,rhi;
 199   rvec       *x;
 200   real       **mdmat,*resnr,**totmdmat;
 201   int        **nmat,**totnmat;
 202   real       *mean_n;
 203   int        *tot_n;
 204   matrix     box;
 205   output_env_t oenv;
 206   gmx_rmpbc_t  gpbc=NULL;
 207
 208   CopyRight(stderr,argv[0]);
 209
 210   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,
 211                     asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 212
 213   fprintf(stderr,"Will truncate at %f nm\n",truncate);
 214   bCalcN = opt2bSet("-no",NFILE,fnm);
 215   bFrames= opt2bSet("-frames",NFILE,fnm);
 216   if ( bCalcN )
 217     fprintf(stderr,"Will calculate number of different contacts\n");
 218
 219   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,FALSE);
 220
 221   fprintf(stderr,"Select group for analysis\n");
 222   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 223
 224   natoms=isize;
 225   snew(useatoms.atom,natoms);
 226   snew(useatoms.atomname,natoms);
 227
 228   useatoms.nres = 0;
 229   snew(useatoms.resinfo,natoms);
 230
 231   prevres = top.atoms.atom[index[0]].resind;
 232   newres  = 0;
 233   for(i=0;(i<isize);i++) {
 234     int ii = index[i];
 235     useatoms.atomname[i]=top.atoms.atomname[ii];
 236     if (top.atoms.atom[ii].resind != prevres) {
 237       prevres = top.atoms.atom[ii].resind;
 238       newres++;
 239       useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 240       if (debug) {
 241         fprintf(debug,"New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 242                 *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 243                 *(top.atoms.atomname[ii]),
 244                 ii,i,newres);
 245       }
 246     }
 247     useatoms.atom[i].resind = newres;
 248   }
 249   useatoms.nres = newres+1;
 250   useatoms.nr = isize;
 251
 252   rndx=res_ndx(&(useatoms));
 253   natm=res_natm(&(useatoms));
 254   nres=useatoms.nres;
 255   fprintf(stderr,"There are %d residues with %d atoms\n",nres,natoms);
 256
 257   snew(resnr,nres);
 258   snew(mdmat,nres);
 259   snew(nmat,nres);
 260   snew(totnmat,nres);
 261   snew(mean_n,nres);
 262   snew(tot_n,nres);
 263   for(i=0; (i<nres); i++) {
 264     snew(mdmat[i],nres);
 265     snew(nmat[i],natoms);
 266     snew(totnmat[i],natoms);
 267     resnr[i]=i+1;
 268   }
 269   snew(totmdmat,nres);
 270   for(i=0; (i<nres); i++)
 271     snew(totmdmat[i],nres);
 272
 273   trxnat=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 274
 275   nframes=0;
 276
 277   rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
 278   rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
 279
 280   gpbc = gmx_rmpbc_init(&top.idef,ePBC,trxnat,box);
 281
 282   if (bFrames)
 283     out=opt2FILE("-frames",NFILE,fnm,"w");
 284   do {
 285     gmx_rmpbc(gpbc,trxnat,box,x);
 286     nframes++;
 287     calc_mat(nres,natoms,rndx,x,index,truncate,mdmat,nmat,ePBC,box);
 288     for (i=0; (i<nres); i++)
 289       for (j=0; (j<natoms); j++)
 290         if (nmat[i][j])
 291           totnmat[i][j]++;
 292     for (i=0; (i<nres); i++)
 293       for (j=0; (j<nres); j++)
 294         totmdmat[i][j] += mdmat[i][j];
 295     if (bFrames) {
 296       sprintf(label,"t=%.0f ps",t);
 297       write_xpm(out,0,label,"Distance (nm)","Residue Index","Residue Index",
 298                 nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 299     }
 300   } while (read_next_x(oenv,status,&t,trxnat,x,box));
 301   fprintf(stderr,"\n");
 302   close_trj(status);
 303   gmx_rmpbc_done(gpbc);
 304   if (bFrames)
 305     ffclose(out);
 306
 307   fprintf(stderr,"Processed %d frames\n",nframes);
 308
 309   for (i=0; (i<nres); i++)
 310     for (j=0; (j<nres); j++)
 311       totmdmat[i][j] /= nframes;
 312   write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),0,"Mean smallest distance",
 313             "Distance (nm)","Residue Index","Residue Index",
 314             nres,nres,resnr,resnr,totmdmat,0,truncate,rlo,rhi,&nlevels);
 315
 316   if ( bCalcN ) {
 317     tot_nmat(nres,natoms,nframes,totnmat,tot_n,mean_n);
 318     fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 319                 "Increase in number of contacts","Residue","Ratio",oenv);
 320     fprintf(fp,"@ legend on\n");
 321     fprintf(fp,"@ legend box on\n");
 322     fprintf(fp,"@ legend loctype view\n");
 323     fprintf(fp,"@ legend 0.75, 0.8\n");
 324     fprintf(fp,"@ legend string 0 \"Total/mean\"\n");
 325     fprintf(fp,"@ legend string 1 \"Total\"\n");
 326     fprintf(fp,"@ legend string 2 \"Mean\"\n");
 327     fprintf(fp,"@ legend string 3 \"# atoms\"\n");
 328     fprintf(fp,"@ legend string 4 \"Mean/# atoms\"\n");
 329     fprintf(fp,"#%3s %8s  %3s  %8s  %3s  %8s\n",
 330             "res","ratio","tot","mean","natm","mean/atm");
 331     for (i=0; (i<nres); i++) {
 332       if (mean_n[i]==0)
 333         ratio=1;
 334       else
 335         ratio=tot_n[i]/mean_n[i];
 336       fprintf(fp,"%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 337               i+1,ratio,tot_n[i],mean_n[i],natm[i],mean_n[i]/natm[i]);
 338     }
 339     ffclose(fp);
 340   }
 341
 342   thanx(stderr);
 343
 344   return 0;
 345 }