src/gromacs/gmxana/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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "macros.h"
  45 #include "vec.h"
  46 #include "sysstuff.h"
  47 #include "typedefs.h"
  48 #include "gromacs/fileio/filenm.h"
  49 #include "gromacs/commandline/pargs.h"
  50 #include "gromacs/fileio/futil.h"
  51 #include "gmx_fatal.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "gromacs/fileio/matio.h"
  54 #include "xvgr.h"
  55 #include "index.h"
  56 #include "gromacs/fileio/tpxio.h"
  57 #include "gromacs/fileio/trxio.h"
  58 #include "rmpbc.h"
  59 #include "pbc.h"
  60 #include "gmx_ana.h"
  61
  62
  63 #define FARAWAY 10000
  64
  65 int *res_ndx(t_atoms *atoms)
  66 {
  67     int *rndx;
  68     int  i, r0;
  69
  70     if (atoms->nr <= 0)
  71     {
  72         return NULL;
  73     }
  74     snew(rndx, atoms->nr);
  75     r0 = atoms->atom[0].resind;
  76     for (i = 0; (i < atoms->nr); i++)
  77     {
  78         rndx[i] = atoms->atom[i].resind-r0;
  79     }
  80
  81     return rndx;
  82 }
  83
  84 int *res_natm(t_atoms *atoms)
  85 {
  86     int *natm;
  87     int  i, j, r0;
  88
  89     if (atoms->nr <= 0)
  90     {
  91         return NULL;
  92     }
  93     snew(natm, atoms->nres);
  94     r0 = atoms->atom[0].resind;
  95     j  = 0;
  96     for (i = 0; (i < atoms->nres); i++)
  97     {
  98         while ((atoms->atom[j].resind)-r0 == i)
  99         {
 100             natm[i]++;
 101             j++;
 102         }
 103     }
 104
 105     return natm;
 106 }
 107
 108 static void calc_mat(int nres, int natoms, int rndx[],
 109                      rvec x[], atom_id *index,
 110                      real trunc, real **mdmat, int **nmat, int ePBC, matrix box)
 111 {
 112     int   i, j, resi, resj;
 113     real  trunc2, r, r2;
 114     t_pbc pbc;
 115     rvec  ddx;
 116
 117     set_pbc(&pbc, ePBC, box);
 118     trunc2 = sqr(trunc);
 119     for (resi = 0; (resi < nres); resi++)
 120     {
 121         for (resj = 0; (resj < nres); resj++)
 122         {
 123             mdmat[resi][resj] = FARAWAY;
 124         }
 125     }
 126     for (i = 0; (i < natoms); i++)
 127     {
 128         resi = rndx[i];
 129         for (j = i+1; (j < natoms); j++)
 130         {
 131             resj = rndx[j];
 132             pbc_dx(&pbc, x[index[i]], x[index[j]], ddx);
 133             r2 = norm2(ddx);
 134             if (r2 < trunc2)
 135             {
 136                 nmat[resi][j]++;
 137                 nmat[resj][i]++;
 138             }
 139             mdmat[resi][resj] = min(r2, mdmat[resi][resj]);
 140         }
 141     }
 142
 143     for (resi = 0; (resi < nres); resi++)
 144     {
 145         mdmat[resi][resi] = 0;
 146         for (resj = resi+1; (resj < nres); resj++)
 147         {
 148             r                 = sqrt(mdmat[resi][resj]);
 149             mdmat[resi][resj] = r;
 150             mdmat[resj][resi] = r;
 151         }
 152     }
 153 }
 154
 155 static void tot_nmat(int nres, int natoms, int nframes, int **nmat,
 156                      int *tot_n, real *mean_n)
 157 {
 158     int i, j;
 159
 160     for (i = 0; (i < nres); i++)
 161     {
 162         for (j = 0; (j < natoms); j++)
 163         {
 164             if (nmat[i][j] != 0)
 165             {
 166                 tot_n[i]++;
 167                 mean_n[i] += nmat[i][j];
 168             }
 169         }
 170         mean_n[i] /= nframes;
 171     }
 172 }
 173
 174 int gmx_mdmat(int argc, char *argv[])
 175 {
 176     const char     *desc[] = {
 177         "[THISMODULE] makes distance matrices consisting of the smallest distance",
 178         "between residue pairs. With [TT]-frames[tt], these distance matrices can be",
 179         "stored in order to see differences in tertiary structure as a",
 180         "function of time. If you choose your options unwisely, this may generate",
 181         "a large output file. By default, only an averaged matrix over the whole",
 182         "trajectory is output.",
 183         "Also a count of the number of different atomic contacts between",
 184         "residues over the whole trajectory can be made.",
 185         "The output can be processed with [gmx-xpm2ps] to make a PostScript (tm) plot."
 186     };
 187     static real     truncate = 1.5;
 188     static gmx_bool bAtom    = FALSE;
 189     static int      nlevels  = 40;
 190     t_pargs         pa[]     = {
 191         { "-t",   FALSE, etREAL, {&truncate},
 192           "trunc distance" },
 193         { "-nlevels",   FALSE, etINT,  {&nlevels},
 194           "Discretize distance in this number of levels" }
 195     };
 196     t_filenm        fnm[] = {
 197         { efTRX, "-f",  NULL, ffREAD },
 198         { efTPS, NULL,  NULL, ffREAD },
 199         { efNDX, NULL,  NULL, ffOPTRD },
 200         { efXPM, "-mean", "dm", ffWRITE },
 201         { efXPM, "-frames", "dmf", ffOPTWR },
 202         { efXVG, "-no", "num", ffOPTWR },
 203     };
 204 #define NFILE asize(fnm)
 205
 206     FILE          *out = NULL, *fp;
 207     t_topology     top;
 208     int            ePBC;
 209     t_atoms        useatoms;
 210     int            isize;
 211     atom_id       *index;
 212     char          *grpname;
 213     int           *rndx, *natm, prevres, newres;
 214
 215     int            i, j, nres, natoms, nframes, it, trxnat;
 216     t_trxstatus   *status;
 217     int            nr0;
 218     gmx_bool       bCalcN, bFrames;
 219     real           t, ratio;
 220     char           title[256], label[234];
 221     t_rgb          rlo, rhi;
 222     rvec          *x;
 223     real         **mdmat, *resnr, **totmdmat;
 224     int          **nmat, **totnmat;
 225     real          *mean_n;
 226     int           *tot_n;
 227     matrix         box;
 228     output_env_t   oenv;
 229     gmx_rmpbc_t    gpbc = NULL;
 230
 231     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE, NFILE, fnm,
 232                            asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 233     {
 234         return 0;
 235     }
 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);
 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, x, box));
 348     fprintf(stderr, "\n");
 349     close_trj(status);
 350     gmx_rmpbc_done(gpbc);
 351     if (bFrames)
 352     {
 353         gmx_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         fprintf(fp, "@ legend on\n");
 375         fprintf(fp, "@ legend box on\n");
 376         fprintf(fp, "@ legend loctype view\n");
 377         fprintf(fp, "@ legend 0.75, 0.8\n");
 378         fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
 379         fprintf(fp, "@ legend string 1 \"Total\"\n");
 380         fprintf(fp, "@ legend string 2 \"Mean\"\n");
 381         fprintf(fp, "@ legend string 3 \"# atoms\"\n");
 382         fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
 383         fprintf(fp, "#%3s %8s  %3s  %8s  %3s  %8s\n",
 384                 "res", "ratio", "tot", "mean", "natm", "mean/atm");
 385         for (i = 0; (i < nres); i++)
 386         {
 387             if (mean_n[i] == 0)
 388             {
 389                 ratio = 1;
 390             }
 391             else
 392             {
 393                 ratio = tot_n[i]/mean_n[i];
 394             }
 395             fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 396                     i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);
 397         }
 398         gmx_ffclose(fp);
 399     }
 400
 401     return 0;
 402 }