src/gromacs/gmxana/gmx_mdmat.cpp

   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 <cmath>
  40 #include <cstring>
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/commandline/filenm.h"
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/fileio/confio.h"
  47 #include "gromacs/fileio/matio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/math/functions.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/pbcutil/rmpbc.h"
  55 #include "gromacs/topology/index.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/utility/arraysize.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/futil.h"
  61 #include "gromacs/utility/smalloc.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[], int *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 = gmx::square(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] = std::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                 = std::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         "[THISMODULE] 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 [gmx-xpm2ps] to make a PostScript (tm) plot."
 187     };
 188     static real     truncate = 1.5;
 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     int              *index;
 212     char             *grpname;
 213     int              *rndx, *natm, prevres, newres;
 214
 215     int               i, j, nres, natoms, nframes, trxnat;
 216     t_trxstatus      *status;
 217     gmx_bool          bCalcN, bFrames;
 218     real              t, ratio;
 219     char              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     gmx_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), &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         xvgrclose(fp);
 400     }
 401
 402     return 0;
 403 }