src/tools/gmx_mdmat.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <string.h>
  44
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "sysstuff.h"
  48 #include "typedefs.h"
  49 #include "filenm.h"
  50 #include "statutil.h"
  51 #include "copyrite.h"
  52 #include "futil.h"
  53 #include "gmx_fatal.h"
  54 #include "smalloc.h"
  55 #include "matio.h"
  56 #include "xvgr.h"
  57 #include "index.h"
  58 #include "tpxio.h"
  59 #include "rmpbc.h"
  60 #include "pbc.h"
  61 #include "gmx_ana.h"
  62
  63
  64 #define FARAWAY 10000
  65
  66 int *res_ndx(t_atoms *atoms)
  67 {
  68   int *rndx;
  69   int i,r0;
  70
  71   if (atoms->nr <= 0)
  72     return NULL;
  73   snew(rndx,atoms->nr);
  74   r0=atoms->atom[0].resind;
  75   for(i=0; (i<atoms->nr); i++)
  76     rndx[i]=atoms->atom[i].resind-r0;
  77
  78   return rndx;
  79 }
  80
  81 int *res_natm(t_atoms *atoms)
  82 {
  83   int *natm;
  84   int i,j,r0;
  85
  86   if (atoms->nr <= 0)
  87     return NULL;
  88   snew(natm,atoms->nres);
  89   r0=atoms->atom[0].resind;
  90   j=0;
  91   for(i=0; (i<atoms->nres); i++) {
  92     while ((atoms->atom[j].resind)-r0 == i) {
  93       natm[i]++;
  94       j++;
  95     }
  96   }
  97
  98   return natm;
  99 }
 100
 101 static void calc_mat(int nres, int natoms, int rndx[],
 102                      rvec x[], atom_id *index,
 103                      real trunc, real **mdmat,int **nmat,int ePBC,matrix box)
 104 {
 105   int i,j,resi,resj;
 106   real trunc2,r,r2;
 107   t_pbc pbc;
 108   rvec ddx;
 109
 110   set_pbc(&pbc,ePBC,box);
 111   trunc2=sqr(trunc);
 112   for(resi=0; (resi<nres); resi++)
 113     for(resj=0; (resj<nres); resj++)
 114       mdmat[resi][resj]=FARAWAY;
 115   for(i=0; (i<natoms); i++) {
 116     resi=rndx[i];
 117     for(j=i+1; (j<natoms); j++) {
 118       resj=rndx[j];
 119       pbc_dx(&pbc,x[index[i]],x[index[j]],ddx);
 120       r2 = norm2(ddx);
 121       if ( r2 < trunc2 ) {
 122         nmat[resi][j]++;
 123         nmat[resj][i]++;
 124       }
 125       mdmat[resi][resj]=min(r2,mdmat[resi][resj]);
 126     }
 127   }
 128
 129   for(resi=0; (resi<nres); resi++) {
 130     mdmat[resi][resi]=0;
 131     for(resj=resi+1; (resj<nres); resj++) {
 132       r=sqrt(mdmat[resi][resj]);
 133       mdmat[resi][resj]=r;
 134       mdmat[resj][resi]=r;
 135     }
 136   }
 137 }
 138
 139 static void tot_nmat(int nres, int natoms, int nframes, int **nmat,
 140                      int *tot_n, real *mean_n)
 141 {
 142   int i,j;
 143
 144   for (i=0; (i<nres); i++) {
 145     for (j=0; (j<natoms); j++)
 146       if (nmat[i][j] != 0) {
 147         tot_n[i]++;
 148         mean_n[i]+=nmat[i][j];
 149       }
 150     mean_n[i]/=nframes;
 151   }
 152 }
 153
 154 int gmx_mdmat(int argc,char *argv[])
 155 {
 156   const char *desc[] = {
 157     "[TT]g_mdmat[tt] makes distance matrices consisting of the smallest distance",
 158     "between residue pairs. With [TT]-frames[tt], these distance matrices can be",
 159     "stored in order to see differences in tertiary structure as a",
 160     "function of time. If you choose your options unwisely, this may generate",
 161     "a large output file. By default, only an averaged matrix over the whole",
 162     "trajectory is output.",
 163     "Also a count of the number of different atomic contacts between",
 164     "residues over the whole trajectory can be made.",
 165     "The output can be processed with [TT]xpm2ps[tt] to make a PostScript (tm) plot."
 166   };
 167   static real truncate=1.5;
 168   static gmx_bool bAtom=FALSE;
 169   static int  nlevels=40;
 170   t_pargs pa[] = {
 171     { "-t",   FALSE, etREAL, {&truncate},
 172       "trunc distance" },
 173     { "-nlevels",   FALSE, etINT,  {&nlevels},
 174       "Discretize distance in this number of levels" }
 175   };
 176   t_filenm   fnm[] = {
 177     { efTRX, "-f",  NULL, ffREAD },
 178     { efTPS, NULL,  NULL, ffREAD },
 179     { efNDX, NULL,  NULL, ffOPTRD },
 180     { efXPM, "-mean", "dm", ffWRITE },
 181     { efXPM, "-frames", "dmf", ffOPTWR },
 182     { efXVG, "-no", "num",ffOPTWR },
 183   };
 184 #define NFILE asize(fnm)
 185
 186   FILE       *out=NULL,*fp;
 187   t_topology top;
 188   int        ePBC;
 189   t_atoms    useatoms;
 190   int        isize;
 191   atom_id    *index;
 192   char       *grpname;
 193   int        *rndx,*natm,prevres,newres;
 194
 195   int        i,j,nres,natoms,nframes,it,trxnat;
 196   t_trxstatus *status;
 197   int        nr0;
 198   gmx_bool       bCalcN,bFrames;
 199   real       t,ratio;
 200   char       title[256],label[234];
 201   t_rgb      rlo,rhi;
 202   rvec       *x;
 203   real       **mdmat,*resnr,**totmdmat;
 204   int        **nmat,**totnmat;
 205   real       *mean_n;
 206   int        *tot_n;
 207   matrix     box;
 208   output_env_t oenv;
 209   gmx_rmpbc_t  gpbc=NULL;
 210
 211   CopyRight(stderr,argv[0]);
 212
 213   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,NFILE,fnm,
 214                     asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 215
 216   fprintf(stderr,"Will truncate at %f nm\n",truncate);
 217   bCalcN = opt2bSet("-no",NFILE,fnm);
 218   bFrames= opt2bSet("-frames",NFILE,fnm);
 219   if ( bCalcN )
 220     fprintf(stderr,"Will calculate number of different contacts\n");
 221
 222   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&x,NULL,box,FALSE);
 223
 224   fprintf(stderr,"Select group for analysis\n");
 225   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 226
 227   natoms=isize;
 228   snew(useatoms.atom,natoms);
 229   snew(useatoms.atomname,natoms);
 230
 231   useatoms.nres = 0;
 232   snew(useatoms.resinfo,natoms);
 233
 234   prevres = top.atoms.atom[index[0]].resind;
 235   newres  = 0;
 236   for(i=0;(i<isize);i++) {
 237     int ii = index[i];
 238     useatoms.atomname[i]=top.atoms.atomname[ii];
 239     if (top.atoms.atom[ii].resind != prevres) {
 240       prevres = top.atoms.atom[ii].resind;
 241       newres++;
 242       useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 243       if (debug) {
 244         fprintf(debug,"New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 245                 *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 246                 *(top.atoms.atomname[ii]),
 247                 ii,i,newres);
 248       }
 249     }
 250     useatoms.atom[i].resind = newres;
 251   }
 252   useatoms.nres = newres+1;
 253   useatoms.nr = isize;
 254
 255   rndx=res_ndx(&(useatoms));
 256   natm=res_natm(&(useatoms));
 257   nres=useatoms.nres;
 258   fprintf(stderr,"There are %d residues with %d atoms\n",nres,natoms);
 259
 260   snew(resnr,nres);
 261   snew(mdmat,nres);
 262   snew(nmat,nres);
 263   snew(totnmat,nres);
 264   snew(mean_n,nres);
 265   snew(tot_n,nres);
 266   for(i=0; (i<nres); i++) {
 267     snew(mdmat[i],nres);
 268     snew(nmat[i],natoms);
 269     snew(totnmat[i],natoms);
 270     resnr[i]=i+1;
 271   }
 272   snew(totmdmat,nres);
 273   for(i=0; (i<nres); i++)
 274     snew(totmdmat[i],nres);
 275
 276   trxnat=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 277
 278   nframes=0;
 279
 280   rlo.r=1.0, rlo.g=1.0, rlo.b=1.0;
 281   rhi.r=0.0, rhi.g=0.0, rhi.b=0.0;
 282
 283   gpbc = gmx_rmpbc_init(&top.idef,ePBC,trxnat,box);
 284
 285   if (bFrames)
 286     out=opt2FILE("-frames",NFILE,fnm,"w");
 287   do {
 288     gmx_rmpbc(gpbc,trxnat,box,x);
 289     nframes++;
 290     calc_mat(nres,natoms,rndx,x,index,truncate,mdmat,nmat,ePBC,box);
 291     for (i=0; (i<nres); i++)
 292       for (j=0; (j<natoms); j++)
 293         if (nmat[i][j])
 294           totnmat[i][j]++;
 295     for (i=0; (i<nres); i++)
 296       for (j=0; (j<nres); j++)
 297         totmdmat[i][j] += mdmat[i][j];
 298     if (bFrames) {
 299       sprintf(label,"t=%.0f ps",t);
 300       write_xpm(out,0,label,"Distance (nm)","Residue Index","Residue Index",
 301                 nres,nres,resnr,resnr,mdmat,0,truncate,rlo,rhi,&nlevels);
 302     }
 303   } while (read_next_x(oenv,status,&t,trxnat,x,box));
 304   fprintf(stderr,"\n");
 305   close_trj(status);
 306   gmx_rmpbc_done(gpbc);
 307   if (bFrames)
 308     ffclose(out);
 309
 310   fprintf(stderr,"Processed %d frames\n",nframes);
 311
 312   for (i=0; (i<nres); i++)
 313     for (j=0; (j<nres); j++)
 314       totmdmat[i][j] /= nframes;
 315   write_xpm(opt2FILE("-mean",NFILE,fnm,"w"),0,"Mean smallest distance",
 316             "Distance (nm)","Residue Index","Residue Index",
 317             nres,nres,resnr,resnr,totmdmat,0,truncate,rlo,rhi,&nlevels);
 318
 319   if ( bCalcN ) {
 320     tot_nmat(nres,natoms,nframes,totnmat,tot_n,mean_n);
 321     fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),
 322                 "Increase in number of contacts","Residue","Ratio",oenv);
 323     fprintf(fp,"@ legend on\n");
 324     fprintf(fp,"@ legend box on\n");
 325     fprintf(fp,"@ legend loctype view\n");
 326     fprintf(fp,"@ legend 0.75, 0.8\n");
 327     fprintf(fp,"@ legend string 0 \"Total/mean\"\n");
 328     fprintf(fp,"@ legend string 1 \"Total\"\n");
 329     fprintf(fp,"@ legend string 2 \"Mean\"\n");
 330     fprintf(fp,"@ legend string 3 \"# atoms\"\n");
 331     fprintf(fp,"@ legend string 4 \"Mean/# atoms\"\n");
 332     fprintf(fp,"#%3s %8s  %3s  %8s  %3s  %8s\n",
 333             "res","ratio","tot","mean","natm","mean/atm");
 334     for (i=0; (i<nres); i++) {
 335       if (mean_n[i]==0)
 336         ratio=1;
 337       else
 338         ratio=tot_n[i]/mean_n[i];
 339       fprintf(fp,"%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 340               i+1,ratio,tot_n[i],mean_n[i],natm[i],mean_n[i]/natm[i]);
 341     }
 342     ffclose(fp);
 343   }
 344
 345   thanx(stderr);
 346
 347   return 0;
 348 }