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,status,nres,natoms,nframes,it,trxnat;
 194   int        nr0;
 195   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
 207   CopyRight(stderr,argv[0]);
 208
 209   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,
 210                     asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 211
 212   fprintf(stderr,"Will truncate at %f nm\n",truncate);
 213   bCalcN = opt2bSet("-no",NFILE,fnm);
 214   bFrames= opt2bSet("-frames",NFILE,fnm);
 215   if ( bCalcN )
 216     fprintf(stderr,"Will calculate number of different contacts\n");
 217
 218   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,FALSE);
 219
 220   fprintf(stderr,"Select group for analysis\n");
 221   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 222
 223   natoms=isize;
 224   snew(useatoms.atom,natoms);
 225   snew(useatoms.atomname,natoms);
 226
 227   useatoms.nres = 0;
 228   snew(useatoms.resinfo,natoms);
 229
 230   prevres = top.atoms.atom[index[0]].resind;
 231   newres  = 0;
 232   for(i=0;(i<isize);i++) {
 233     int ii = index[i];
 234     useatoms.atomname[i]=top.atoms.atomname[ii];
 235     if (top.atoms.atom[ii].resind != prevres) {
 236       prevres = top.atoms.atom[ii].resind;
 237       newres++;
 238       useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 239       if (debug) {
 240         fprintf(debug,"New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 241                 *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 242                 *(top.atoms.atomname[ii]),
 243                 ii,i,newres);
 244       }
 245     }
 246     useatoms.atom[i].resind = newres;
 247   }
 248   useatoms.nres = newres+1;
 249   useatoms.nr = isize;
 250
 251   rndx=res_ndx(&(useatoms));
 252   natm=res_natm(&(useatoms));
 253   nres=useatoms.nres;
 254   fprintf(stderr,"There are %d residues with %d atoms\n",nres,natoms);
 255
 256   snew(resnr,nres);
 257   snew(mdmat,nres);
 258   snew(nmat,nres);
 259   snew(totnmat,nres);
 260   snew(mean_n,nres);
 261   snew(tot_n,nres);
 262   for(i=0; (i<nres); i++) {
 263     snew(mdmat[i],nres);
 264     snew(nmat[i],natoms);
 265     snew(totnmat[i],natoms);
 266     resnr[i]=i+1;
 267   }
 268   snew(totmdmat,nres);
 269   for(i=0; (i<nres); i++)
 270     snew(totmdmat[i],nres);
 271
 272   trxnat=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 273
 274   nframes=0;
 275
 276   rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
 277   rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
 278   if (bFrames)
 279     out=opt2FILE("-frames",NFILE,fnm,"w");
 280   do {
 281     rm_pbc(&top.idef,ePBC,trxnat,box,x,x);
 282     nframes++;
 283     calc_mat(nres,natoms,rndx,x,index,truncate,mdmat,nmat,ePBC,box);
 284     for (i=0; (i<nres); i++)
 285       for (j=0; (j<natoms); j++)
 286         if (nmat[i][j])
 287           totnmat[i][j]++;
 288     for (i=0; (i<nres); i++)
 289       for (j=0; (j<nres); j++)
 290         totmdmat[i][j] += mdmat[i][j];
 291     if (bFrames) {
 292       sprintf(label,"t=%.0f ps",t);
 293       write_xpm(out,0,label,"Distance (nm)","Residue Index","Residue Index",
 294                 nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 295     }
 296   } while (read_next_x(oenv,status,&t,trxnat,x,box));
 297   fprintf(stderr,"\n");
 298   close_trj(status);
 299   if (bFrames)
 300     ffclose(out);
 301
 302   fprintf(stderr,"Processed %d frames\n",nframes);
 303
 304   for (i=0; (i<nres); i++)
 305     for (j=0; (j<nres); j++)
 306       totmdmat[i][j] /= nframes;
 307   write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),0,"Mean smallest distance",
 308             "Distance (nm)","Residue Index","Residue Index",
 309             nres,nres,resnr,resnr,totmdmat,0,truncate,rlo,rhi,&nlevels);
 310
 311   if ( bCalcN ) {
 312     tot_nmat(nres,natoms,nframes,totnmat,tot_n,mean_n);
 313     fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 314                 "Increase in number of contacts","Residue","Ratio",oenv);
 315     fprintf(fp,"@ legend on\n");
 316     fprintf(fp,"@ legend box on\n");
 317     fprintf(fp,"@ legend loctype view\n");
 318     fprintf(fp,"@ legend 0.75, 0.8\n");
 319     fprintf(fp,"@ legend string 0 \"Total/mean\"\n");
 320     fprintf(fp,"@ legend string 1 \"Total\"\n");
 321     fprintf(fp,"@ legend string 2 \"Mean\"\n");
 322     fprintf(fp,"@ legend string 3 \"# atoms\"\n");
 323     fprintf(fp,"@ legend string 4 \"Mean/# atoms\"\n");
 324     fprintf(fp,"#%3s %8s  %3s  %8s  %3s  %8s\n",
 325             "res","ratio","tot","mean","natm","mean/atm");
 326     for (i=0; (i<nres); i++) {
 327       if (mean_n[i]==0)
 328         ratio=1;
 329       else
 330         ratio=tot_n[i]/mean_n[i];
 331       fprintf(fp,"%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 332               i+1,ratio,tot_n[i],mean_n[i],natm[i],mean_n[i]/natm[i]);
 333     }
 334     ffclose(fp);
 335   }
 336
 337   thanx(stderr);
 338
 339   return 0;
 340 }