src/gromacs/gmxana/gmx_mdmat.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41
  42 #include "macros.h"
  43 #include "vec.h"
  44 #include "sysstuff.h"
  45 #include "typedefs.h"
  46 #include "gromacs/fileio/filenm.h"
  47 #include "gromacs/commandline/pargs.h"
  48 #include "gromacs/fileio/futil.h"
  49 #include "gmx_fatal.h"
  50 #include "smalloc.h"
  51 #include "gromacs/fileio/matio.h"
  52 #include "xvgr.h"
  53 #include "index.h"
  54 #include "gromacs/fileio/tpxio.h"
  55 #include "gromacs/fileio/trxio.h"
  56 #include "rmpbc.h"
  57 #include "pbc.h"
  58 #include "gmx_ana.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;
 226     output_env_t   oenv;
 227     gmx_rmpbc_t    gpbc = NULL;
 228
 229     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE, 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         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         tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
 370         fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
 371                       "Increase in number of contacts", "Residue", "Ratio", oenv);
 372         fprintf(fp, "@ legend on\n");
 373         fprintf(fp, "@ legend box on\n");
 374         fprintf(fp, "@ legend loctype view\n");
 375         fprintf(fp, "@ legend 0.75, 0.8\n");
 376         fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
 377         fprintf(fp, "@ legend string 1 \"Total\"\n");
 378         fprintf(fp, "@ legend string 2 \"Mean\"\n");
 379         fprintf(fp, "@ legend string 3 \"# atoms\"\n");
 380         fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
 381         fprintf(fp, "#%3s %8s  %3s  %8s  %3s  %8s\n",
 382                 "res", "ratio", "tot", "mean", "natm", "mean/atm");
 383         for (i = 0; (i < nres); i++)
 384         {
 385             if (mean_n[i] == 0)
 386             {
 387                 ratio = 1;
 388             }
 389             else
 390             {
 391                 ratio = tot_n[i]/mean_n[i];
 392             }
 393             fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 394                     i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);
 395         }
 396         ffclose(fp);
 397     }
 398
 399     return 0;
 400 }