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,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/filenm.h"
  46 #include "gromacs/commandline/pargs.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/matio.h"
  49 #include "gromacs/fileio/trxio.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/pbcutil/pbc.h"
  55 #include "gromacs/pbcutil/rmpbc.h"
  56 #include "gromacs/topology/index.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/arraysize.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64
  65 #define FARAWAY 10000
  66
  67 static int* res_ndx(t_atoms* atoms)
  68 {
  69     int* rndx;
  70     int  i, r0;
  71
  72     if (atoms->nr <= 0)
  73     {
  74         return nullptr;
  75     }
  76     snew(rndx, atoms->nr);
  77     r0 = atoms->atom[0].resind;
  78     for (i = 0; (i < atoms->nr); i++)
  79     {
  80         rndx[i] = atoms->atom[i].resind - r0;
  81     }
  82
  83     return rndx;
  84 }
  85
  86 static int* res_natm(t_atoms* atoms)
  87 {
  88     int* natm;
  89     int  i, j, r0;
  90
  91     if (atoms->nr <= 0)
  92     {
  93         return nullptr;
  94     }
  95     snew(natm, atoms->nres);
  96     r0 = atoms->atom[0].resind;
  97     j  = 0;
  98     for (i = 0; (i < atoms->nres); i++)
  99     {
 100         while ((atoms->atom[j].resind) - r0 == i)
 101         {
 102             natm[i]++;
 103             j++;
 104         }
 105     }
 106
 107     return natm;
 108 }
 109
 110 static void calc_mat(int        nres,
 111                      int        natoms,
 112                      const int  rndx[],
 113                      rvec       x[],
 114                      const int* index,
 115                      real       trunc,
 116                      real**     mdmat,
 117                      int**      nmat,
 118                      int        ePBC,
 119                      matrix     box)
 120 {
 121     int   i, j, resi, resj;
 122     real  trunc2, r, r2;
 123     t_pbc pbc;
 124     rvec  ddx;
 125
 126     set_pbc(&pbc, ePBC, box);
 127     trunc2 = gmx::square(trunc);
 128     for (resi = 0; (resi < nres); resi++)
 129     {
 130         for (resj = 0; (resj < nres); resj++)
 131         {
 132             mdmat[resi][resj] = FARAWAY;
 133         }
 134     }
 135     for (i = 0; (i < natoms); i++)
 136     {
 137         resi = rndx[i];
 138         for (j = i + 1; (j < natoms); j++)
 139         {
 140             resj = rndx[j];
 141             pbc_dx(&pbc, x[index[i]], x[index[j]], ddx);
 142             r2 = norm2(ddx);
 143             if (r2 < trunc2)
 144             {
 145                 nmat[resi][j]++;
 146                 nmat[resj][i]++;
 147             }
 148             mdmat[resi][resj] = std::min(r2, mdmat[resi][resj]);
 149         }
 150     }
 151
 152     for (resi = 0; (resi < nres); resi++)
 153     {
 154         mdmat[resi][resi] = 0;
 155         for (resj = resi + 1; (resj < nres); resj++)
 156         {
 157             r                 = std::sqrt(mdmat[resi][resj]);
 158             mdmat[resi][resj] = r;
 159             mdmat[resj][resi] = r;
 160         }
 161     }
 162 }
 163
 164 static void tot_nmat(int nres, int natoms, int nframes, int** nmat, int* tot_n, real* mean_n)
 165 {
 166     int i, j;
 167
 168     for (i = 0; (i < nres); i++)
 169     {
 170         for (j = 0; (j < natoms); j++)
 171         {
 172             if (nmat[i][j] != 0)
 173             {
 174                 tot_n[i]++;
 175                 mean_n[i] += nmat[i][j];
 176             }
 177         }
 178         mean_n[i] /= nframes;
 179     }
 180 }
 181
 182 int gmx_mdmat(int argc, char* argv[])
 183 {
 184     const char* desc[] = {
 185         "[THISMODULE] makes distance matrices consisting of the smallest distance",
 186         "between residue pairs. With [TT]-frames[tt], these distance matrices can be",
 187         "stored in order to see differences in tertiary structure as a",
 188         "function of time. If you choose your options unwisely, this may generate",
 189         "a large output file. By default, only an averaged matrix over the whole",
 190         "trajectory is output.",
 191         "Also a count of the number of different atomic contacts between",
 192         "residues over the whole trajectory can be made.",
 193         "The output can be processed with [gmx-xpm2ps] to make a PostScript (tm) plot."
 194     };
 195     static real truncate = 1.5;
 196     static int  nlevels  = 40;
 197     t_pargs     pa[]     = {
 198         { "-t", FALSE, etREAL, { &truncate }, "trunc distance" },
 199         { "-nlevels", FALSE, etINT, { &nlevels }, "Discretize distance in this number of levels" }
 200     };
 201     t_filenm fnm[] = {
 202         { efTRX, "-f", nullptr, ffREAD },     { efTPS, nullptr, nullptr, ffREAD },
 203         { efNDX, nullptr, nullptr, ffOPTRD }, { efXPM, "-mean", "dm", ffWRITE },
 204         { efXPM, "-frames", "dmf", ffOPTWR }, { efXVG, "-no", "num", ffOPTWR },
 205     };
 206 #define NFILE asize(fnm)
 207
 208     FILE *     out = nullptr, *fp;
 209     t_topology top;
 210     int        ePBC;
 211     t_atoms    useatoms;
 212     int        isize;
 213     int*       index;
 214     char*      grpname;
 215     int *      rndx, *natm, prevres, newres;
 216
 217     int               i, j, nres, natoms, nframes, trxnat;
 218     t_trxstatus*      status;
 219     gmx_bool          bCalcN, bFrames;
 220     real              t, ratio;
 221     char              label[234];
 222     t_rgb             rlo, rhi;
 223     rvec*             x;
 224     real **           mdmat, *resnr, **totmdmat;
 225     int **            nmat, **totnmat;
 226     real*             mean_n;
 227     int*              tot_n;
 228     matrix            box = { { 0 } };
 229     gmx_output_env_t* oenv;
 230     gmx_rmpbc_t       gpbc = nullptr;
 231
 232     if (!parse_common_args(&argc, argv, PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 233                            0, nullptr, &oenv))
 234     {
 235         return 0;
 236     }
 237
 238     fprintf(stderr, "Will truncate at %f nm\n", truncate);
 239     bCalcN  = opt2bSet("-no", NFILE, fnm);
 240     bFrames = opt2bSet("-frames", NFILE, fnm);
 241     if (bCalcN)
 242     {
 243         fprintf(stderr, "Will calculate number of different contacts\n");
 244     }
 245
 246     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &x, nullptr, box, FALSE);
 247
 248     fprintf(stderr, "Select group for analysis\n");
 249     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 250
 251     natoms = isize;
 252     snew(useatoms.atom, natoms);
 253     snew(useatoms.atomname, natoms);
 254
 255     useatoms.nres = 0;
 256     snew(useatoms.resinfo, natoms);
 257
 258     prevres = top.atoms.atom[index[0]].resind;
 259     newres  = 0;
 260     for (i = 0; (i < isize); i++)
 261     {
 262         int ii               = index[i];
 263         useatoms.atomname[i] = top.atoms.atomname[ii];
 264         if (top.atoms.atom[ii].resind != prevres)
 265         {
 266             prevres = top.atoms.atom[ii].resind;
 267             newres++;
 268             useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 269             if (debug)
 270             {
 271                 fprintf(debug, "New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 272                         *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 273                         *(top.atoms.atomname[ii]), 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;
 310     rlo.g = 1.0;
 311     rlo.b = 1.0;
 312     rhi.r = 0.0;
 313     rhi.g = 0.0;
 314     rhi.b = 0.0;
 315
 316     gpbc = gmx_rmpbc_init(&top.idef, ePBC, trxnat);
 317
 318     if (bFrames)
 319     {
 320         out = opt2FILE("-frames", NFILE, fnm, "w");
 321     }
 322     do
 323     {
 324         gmx_rmpbc(gpbc, trxnat, box, x);
 325         nframes++;
 326         calc_mat(nres, natoms, rndx, x, index, truncate, mdmat, nmat, ePBC, box);
 327         for (i = 0; (i < nres); i++)
 328         {
 329             for (j = 0; (j < natoms); j++)
 330             {
 331                 if (nmat[i][j])
 332                 {
 333                     totnmat[i][j]++;
 334                 }
 335             }
 336         }
 337         for (i = 0; (i < nres); i++)
 338         {
 339             for (j = 0; (j < nres); j++)
 340             {
 341                 totmdmat[i][j] += mdmat[i][j];
 342             }
 343         }
 344         if (bFrames)
 345         {
 346             sprintf(label, "t=%.0f ps", t);
 347             write_xpm(out, 0, label, "Distance (nm)", "Residue Index", "Residue Index", nres, nres,
 348                       resnr, resnr, mdmat, 0, truncate, rlo, rhi, &nlevels);
 349         }
 350     } while (read_next_x(oenv, status, &t, x, box));
 351     fprintf(stderr, "\n");
 352     close_trx(status);
 353     gmx_rmpbc_done(gpbc);
 354     if (bFrames)
 355     {
 356         gmx_ffclose(out);
 357     }
 358
 359     fprintf(stderr, "Processed %d frames\n", nframes);
 360
 361     for (i = 0; (i < nres); i++)
 362     {
 363         for (j = 0; (j < nres); j++)
 364         {
 365             totmdmat[i][j] /= nframes;
 366         }
 367     }
 368     write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean smallest distance", "Distance (nm)",
 369               "Residue Index", "Residue Index", nres, nres, resnr, resnr, totmdmat, 0, truncate,
 370               rlo, rhi, &nlevels);
 371
 372     if (bCalcN)
 373     {
 374         char** legend;
 375
 376         snew(legend, 5);
 377         for (i = 0; i < 5; i++)
 378         {
 379             snew(legend[i], STRLEN);
 380         }
 381         tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
 382         fp = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Increase in number of contacts", "Residue",
 383                       "Ratio", oenv);
 384         sprintf(legend[0], "Total/mean");
 385         sprintf(legend[1], "Total");
 386         sprintf(legend[2], "Mean");
 387         sprintf(legend[3], "# atoms");
 388         sprintf(legend[4], "Mean/# atoms");
 389         xvgr_legend(fp, 5, legend, oenv);
 390         for (i = 0; (i < nres); i++)
 391         {
 392             if (mean_n[i] == 0)
 393             {
 394                 ratio = 1;
 395             }
 396             else
 397             {
 398                 ratio = tot_n[i] / mean_n[i];
 399             }
 400             fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n", i + 1, ratio, tot_n[i], mean_n[i],
 401                     natm[i], mean_n[i] / natm[i]);
 402         }
 403         xvgrclose(fp);
 404     }
 405
 406     return 0;
 407 }