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