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