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