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