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