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