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