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     if (!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         return 0;
 232     }
 233
 234     fprintf(stderr, "Will truncate at %f nm\n", truncate);
 235     bCalcN  = opt2bSet("-no", NFILE, fnm);
 236     bFrames = opt2bSet("-frames", NFILE, fnm);
 237     if (bCalcN)
 238     {
 239         fprintf(stderr, "Will calculate number of different contacts\n");
 240     }
 241
 242     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box, FALSE);
 243
 244     fprintf(stderr, "Select group for analysis\n");
 245     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 246
 247     natoms = isize;
 248     snew(useatoms.atom, natoms);
 249     snew(useatoms.atomname, natoms);
 250
 251     useatoms.nres = 0;
 252     snew(useatoms.resinfo, natoms);
 253
 254     prevres = top.atoms.atom[index[0]].resind;
 255     newres  = 0;
 256     for (i = 0; (i < isize); i++)
 257     {
 258         int ii = index[i];
 259         useatoms.atomname[i] = top.atoms.atomname[ii];
 260         if (top.atoms.atom[ii].resind != prevres)
 261         {
 262             prevres = top.atoms.atom[ii].resind;
 263             newres++;
 264             useatoms.resinfo[i] = top.atoms.resinfo[prevres];
 265             if (debug)
 266             {
 267                 fprintf(debug, "New residue: atom %5s %5s %6d, index entry %5d, newres %5d\n",
 268                         *(top.atoms.resinfo[top.atoms.atom[ii].resind].name),
 269                         *(top.atoms.atomname[ii]),
 270                         ii, i, newres);
 271             }
 272         }
 273         useatoms.atom[i].resind = newres;
 274     }
 275     useatoms.nres = newres+1;
 276     useatoms.nr   = isize;
 277
 278     rndx = res_ndx(&(useatoms));
 279     natm = res_natm(&(useatoms));
 280     nres = useatoms.nres;
 281     fprintf(stderr, "There are %d residues with %d atoms\n", nres, natoms);
 282
 283     snew(resnr, nres);
 284     snew(mdmat, nres);
 285     snew(nmat, nres);
 286     snew(totnmat, nres);
 287     snew(mean_n, nres);
 288     snew(tot_n, nres);
 289     for (i = 0; (i < nres); i++)
 290     {
 291         snew(mdmat[i], nres);
 292         snew(nmat[i], natoms);
 293         snew(totnmat[i], natoms);
 294         resnr[i] = i+1;
 295     }
 296     snew(totmdmat, nres);
 297     for (i = 0; (i < nres); i++)
 298     {
 299         snew(totmdmat[i], nres);
 300     }
 301
 302     trxnat = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 303
 304     nframes = 0;
 305
 306     rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
 307     rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;
 308
 309     gpbc = gmx_rmpbc_init(&top.idef, ePBC, trxnat);
 310
 311     if (bFrames)
 312     {
 313         out = opt2FILE("-frames", NFILE, fnm, "w");
 314     }
 315     do
 316     {
 317         gmx_rmpbc(gpbc, trxnat, box, x);
 318         nframes++;
 319         calc_mat(nres, natoms, rndx, x, index, truncate, mdmat, nmat, ePBC, box);
 320         for (i = 0; (i < nres); i++)
 321         {
 322             for (j = 0; (j < natoms); j++)
 323             {
 324                 if (nmat[i][j])
 325                 {
 326                     totnmat[i][j]++;
 327                 }
 328             }
 329         }
 330         for (i = 0; (i < nres); i++)
 331         {
 332             for (j = 0; (j < nres); j++)
 333             {
 334                 totmdmat[i][j] += mdmat[i][j];
 335             }
 336         }
 337         if (bFrames)
 338         {
 339             sprintf(label, "t=%.0f ps", t);
 340             write_xpm(out, 0, label, "Distance (nm)", "Residue Index", "Residue Index",
 341                       nres, nres, resnr, resnr, mdmat, 0, truncate, rlo, rhi, &nlevels);
 342         }
 343     }
 344     while (read_next_x(oenv, status, &t, x, box));
 345     fprintf(stderr, "\n");
 346     close_trj(status);
 347     gmx_rmpbc_done(gpbc);
 348     if (bFrames)
 349     {
 350         ffclose(out);
 351     }
 352
 353     fprintf(stderr, "Processed %d frames\n", nframes);
 354
 355     for (i = 0; (i < nres); i++)
 356     {
 357         for (j = 0; (j < nres); j++)
 358         {
 359             totmdmat[i][j] /= nframes;
 360         }
 361     }
 362     write_xpm(opt2FILE("-mean", NFILE, fnm, "w"), 0, "Mean smallest distance",
 363               "Distance (nm)", "Residue Index", "Residue Index",
 364               nres, nres, resnr, resnr, totmdmat, 0, truncate, rlo, rhi, &nlevels);
 365
 366     if (bCalcN)
 367     {
 368         tot_nmat(nres, natoms, nframes, totnmat, tot_n, mean_n);
 369         fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
 370                       "Increase in number of contacts", "Residue", "Ratio", oenv);
 371         fprintf(fp, "@ legend on\n");
 372         fprintf(fp, "@ legend box on\n");
 373         fprintf(fp, "@ legend loctype view\n");
 374         fprintf(fp, "@ legend 0.75, 0.8\n");
 375         fprintf(fp, "@ legend string 0 \"Total/mean\"\n");
 376         fprintf(fp, "@ legend string 1 \"Total\"\n");
 377         fprintf(fp, "@ legend string 2 \"Mean\"\n");
 378         fprintf(fp, "@ legend string 3 \"# atoms\"\n");
 379         fprintf(fp, "@ legend string 4 \"Mean/# atoms\"\n");
 380         fprintf(fp, "#%3s %8s  %3s  %8s  %3s  %8s\n",
 381                 "res", "ratio", "tot", "mean", "natm", "mean/atm");
 382         for (i = 0; (i < nres); i++)
 383         {
 384             if (mean_n[i] == 0)
 385             {
 386                 ratio = 1;
 387             }
 388             else
 389             {
 390                 ratio = tot_n[i]/mean_n[i];
 391             }
 392             fprintf(fp, "%3d  %8.3f  %3d  %8.3f  %3d  %8.3f\n",
 393                     i+1, ratio, tot_n[i], mean_n[i], natm[i], mean_n[i]/natm[i]);
 394         }
 395         ffclose(fp);
 396     }
 397
 398     return 0;
 399 }