src/tools/gmx_select.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  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \example gmx_select.c
  32  * \brief Utility/example program for writing out basic data for selections.
  33  */
  34 /*! \file
  35  * \brief Utility program for writing out basic data for selections.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <copyrite.h>
  42 #include <index.h>
  43 #include <macros.h>
  44 #include <smalloc.h>
  45 #include <statutil.h>
  46 #include <xvgr.h>
  47 #include <string2.h>
  48 #include <trajana.h>
  49 #include "gmx_ana.h"
  50 #include "gmx_fatal.h"
  51
  52
  53 typedef struct
  54 {
  55     bool                bDump;
  56     bool                bFracNorm;
  57     const char         *routt;
  58     int                *size;
  59     FILE               *sfp;
  60     FILE               *cfp;
  61     FILE               *ifp;
  62     t_blocka           *block;
  63     char              **gnames;
  64     FILE               *mfp;
  65     gmx_ana_indexmap_t *mmap;
  66 } t_dsdata;
  67
  68 static int
  69 print_data(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  70            int nr, gmx_ana_selection_t *sel[], void *data)
  71 {
  72     t_dsdata           *d = (t_dsdata *)data;
  73     int                 g, i, b, mask;
  74     real                normfac;
  75     char                buf2[100],*buf,*nl;
  76     static int          bFirstFrame=1;
  77
  78     /* Write the sizes of the groups, possibly normalized */
  79     if (d->sfp)
  80     {
  81         fprintf(d->sfp, "%11.3f", fr->time);
  82         for (g = 0; g < nr; ++g)
  83         {
  84             normfac = d->bFracNorm ? 1.0 / sel[g]->cfrac : 1.0;
  85             fprintf(d->sfp, " %8.3f", sel[g]->p.nr * normfac / d->size[g]);
  86         }
  87         fprintf(d->sfp, "\n");
  88     }
  89
  90     /* Write the covered fraction */
  91     if (d->cfp)
  92     {
  93         fprintf(d->cfp, "%11.3f", fr->time);
  94         for (g = 0; g < nr; ++g)
  95         {
  96             fprintf(d->cfp, " %6.4f", sel[g]->cfrac);
  97         }
  98         fprintf(d->cfp, "\n");
  99     }
 100
 101     /* Write the actual indices */
 102     if (d->ifp)
 103     {
 104         if (!d->bDump)
 105         {
 106             fprintf(d->ifp, "%11.3f", fr->time);
 107         }
 108         for (g = 0; g < nr; ++g)
 109         {
 110             if (!d->bDump)
 111             {
 112                 fprintf(d->ifp, " %d", sel[g]->p.nr);
 113             }
 114             for (i = 0; i < sel[g]->p.nr; ++i)
 115             {
 116                 if (sel[g]->p.m.type == INDEX_RES && d->routt[0] == 'n')
 117                 {
 118                     fprintf(d->ifp, " %d", top->atoms.resinfo[sel[g]->p.m.mapid[i]].nr);
 119                 }
 120                 else
 121                 {
 122                     fprintf(d->ifp, " %d", sel[g]->p.m.mapid[i]+1);
 123                 }
 124             }
 125         }
 126         fprintf(d->ifp, "\n");
 127     }
 128
 129     if (d->block)
 130     {
 131         for (g = 0; g < nr; ++g)
 132         {
 133             if (sel[g]->bDynamic || bFirstFrame)
 134             {
 135                 buf = strdup(sel[g]->name);
 136                 while ((nl = strchr(buf, ' ')) != NULL)
 137                 {
 138                     *nl = '_';
 139                 }
 140                 if (sel[g]->bDynamic)
 141                 {
 142                     sprintf(buf2, "_%.3f", fr->time);
 143                     srenew(buf, strlen(buf) + strlen(buf2) + 1);
 144                     strcat(buf, buf2);
 145                 }
 146                 add_grp(d->block, &d->gnames, sel[g]->p.nr, sel[g]->p.m.mapid, buf);
 147                 sfree(buf);
 148             }
 149         }
 150     }
 151
 152     /* Write masks */
 153     if (d->mfp)
 154     {
 155         gmx_ana_indexmap_update(d->mmap, sel[0]->g, TRUE);
 156         if (!d->bDump)
 157         {
 158             fprintf(d->mfp, "%11.3f", fr->time);
 159         }
 160         for (b = 0; b < d->mmap->nr; ++b)
 161         {
 162             mask = (d->mmap->refid[b] == -1 ? 0 : 1);
 163             fprintf(d->mfp, d->bDump ? "%d\n" : " %d", mask);
 164         }
 165         if (!d->bDump)
 166         {
 167             fprintf(d->mfp, "\n");
 168         }
 169     }
 170     bFirstFrame = 0;
 171     return 0;
 172 }
 173
 174 int
 175 gmx_select(int argc, char *argv[])
 176 {
 177     const char *desc[] = {
 178         "g_select writes out basic data about dynamic selections.",
 179         "It can be used for some simple analyses, or the output can",
 180         "be combined with output from other programs and/or external",
 181         "analysis programs to calculate more complex things.",
 182         "Any combination of the output options is possible, but note",
 183         "that [TT]-om[tt] only operates on the first selection.[PAR]",
 184         "With [TT]-os[tt], calculates the number of positions in each",
 185         "selection for each frame. With [TT]-norm[tt], the output is",
 186         "between 0 and 1 and describes the fraction from the maximum",
 187         "number of positions (e.g., for selection 'resname RA and x < 5'",
 188         "the maximum number of positions is the number of atoms in",
 189         "RA residues). With [TT]-cfnorm[tt], the output is divided",
 190         "by the fraction covered by the selection.",
 191         "[TT]-norm[tt] and [TT]-cfnorm[tt] can be specified independently",
 192         "of one another.[PAR]",
 193         "With [TT]-oc[tt], the fraction covered by each selection is",
 194         "written out as a function of time.[PAR]",
 195         "With [TT]-oi[tt], the selected atoms/residues/molecules are",
 196         "written out as a function of time. In the output, the first",
 197         "column contains the frame time, the second contains the number",
 198         "of positions, followed by the atom/residue/molecule numbers.",
 199         "If more than one selection is specified, the size of the second",
 200         "group immediately follows the last number of the first group",
 201         "and so on. With [TT]-dump[tt], the frame time and the number",
 202         "of positions is omitted from the output. In this case, only one",
 203         "selection can be given.[PAR]",
 204         "With [TT]-on[tt], the selected atoms are written as a index file",
 205         "compatible with make_ndx and the analyzing tools. Each selection",
 206         "is written as a selection group and for dynamic selections a",
 207         "group is written for each frame.[PAR]",
 208         "For residue numbers, the output of [TT]-oi[tt] can be controlled",
 209         "with [TT]-resnr[tt]: [TT]number[tt] (default) prints the residue",
 210         "numbers as they appear in the input file, while [TT]index[tt] prints",
 211         "unique numbers assigned to the residues in the order they appear",
 212         "in the input file, starting with 1. The former is more intuitive,",
 213         "but if the input contains multiple residues with the same number,",
 214         "the output can be less useful.[PAR]",
 215         "With [TT]-om[tt], a mask is printed for the first selection",
 216         "as a function of time. Each line in the output corresponds to",
 217         "one frame, and contains either 0/1 for each atom/residue/molecule",
 218         "possibly selected. 1 stands for the atom/residue/molecule being",
 219         "selected for the current frame, 0 for not selected.",
 220         "With [TT]-dump[tt], the frame time is omitted from the output.",
 221     };
 222
 223     bool                bDump     = FALSE;
 224     bool                bFracNorm = FALSE;
 225     bool                bTotNorm  = FALSE;
 226     const char         *routt[] = {NULL, "number", "index", NULL};
 227     t_pargs             pa[] = {
 228         {"-dump",   FALSE, etBOOL, {&bDump},
 229          "Do not print the frame time (-om, -oi) or the index size (-oi)"},
 230         {"-norm",   FALSE, etBOOL, {&bTotNorm},
 231          "Normalize by total number of positions with -os"},
 232         {"-cfnorm", FALSE, etBOOL, {&bFracNorm},
 233          "Normalize by covered fraction with -os"},
 234         {"-resnr",  FALSE, etENUM, {routt},
 235          "Residue number output type"},
 236     };
 237
 238     t_filenm            fnm[] = {
 239         {efXVG, "-os", "size.xvg",  ffOPTWR},
 240         {efXVG, "-oc", "cfrac.xvg", ffOPTWR},
 241         {efDAT, "-oi", "index.dat", ffOPTWR},
 242         {efDAT, "-om", "mask.dat",  ffOPTWR},
 243         {efNDX, "-on", "index.ndx", ffOPTWR},
 244     };
 245 #define NFILE asize(fnm)
 246
 247     gmx_ana_traj_t       *trj;
 248     t_topology           *top;
 249     int                   ngrps;
 250     gmx_ana_selection_t **sel;
 251     char                **grpnames;
 252     t_dsdata              d;
 253     const char            *fnSize, *fnFrac, *fnIndex, *fnNdx, *fnMask;
 254     int                   g;
 255     int                   rc;
 256     output_env_t          oenv;
 257
 258     CopyRight(stderr, argv[0]);
 259     gmx_ana_traj_create(&trj, 0);
 260     gmx_ana_set_nanagrps(trj, -1);
 261     parse_trjana_args(trj, &argc, argv, 0,
 262                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
 263                       &oenv);
 264     gmx_ana_get_nanagrps(trj, &ngrps);
 265     gmx_ana_get_anagrps(trj, &sel);
 266     gmx_ana_init_coverfrac(trj, CFRAC_SOLIDANGLE);
 267
 268     /* Get output file names */
 269     fnSize  = opt2fn_null("-os", NFILE, fnm);
 270     fnFrac  = opt2fn_null("-oc", NFILE, fnm);
 271     fnIndex = opt2fn_null("-oi", NFILE, fnm);
 272     fnNdx   = opt2fn_null("-on", NFILE, fnm);
 273     fnMask  = opt2fn_null("-om", NFILE, fnm);
 274     /* Write out sizes if nothing specified */
 275     if (!fnFrac && !fnIndex && !fnMask && !fnNdx)
 276     {
 277         fnSize = opt2fn("-os", NFILE, fnm);
 278     }
 279
 280     if ( bDump && ngrps > 1)
 281     {
 282         gmx_fatal(FARGS, "Only one index group allowed with -dump");
 283     }
 284     if (fnNdx && sel[0]->p.m.type != INDEX_ATOM)
 285     {
 286         gmx_fatal(FARGS, "Only atom selection allowed with -on");
 287     }
 288     if (fnMask && ngrps > 1)
 289     {
 290         fprintf(stderr, "warning: the mask (-om) will only be written for the first group\n");
 291     }
 292     if (fnMask && !sel[0]->bDynamic)
 293     {
 294         fprintf(stderr, "warning: will not write the mask (-om) for a static selection\n");
 295         fnMask = NULL;
 296     }
 297
 298     /* Initialize reference calculation for masks */
 299     if (fnMask)
 300     {
 301         gmx_ana_get_topology(trj, FALSE, &top, NULL);
 302         snew(d.mmap, 1);
 303         gmx_ana_indexmap_init(d.mmap, sel[0]->g, top, sel[0]->p.m.type);
 304     }
 305
 306     /* Initialize calculation data */
 307     d.bDump     = bDump;
 308     d.bFracNorm = bFracNorm;
 309     d.routt     = routt[0];
 310     snew(d.size,  ngrps);
 311     for (g = 0; g < ngrps; ++g)
 312     {
 313         d.size[g] = bTotNorm ? sel[g]->p.nr : 1;
 314     }
 315
 316     /* Open output files */
 317     d.sfp = d.cfp = d.ifp = d.mfp = NULL;
 318     d.block = NULL;
 319     gmx_ana_get_grpnames(trj, &grpnames);
 320     if (fnSize)
 321     {
 322         d.sfp = xvgropen(fnSize, "Selection size", "Time (ps)", "Number",oenv);
 323         xvgr_selections(d.sfp, trj);
 324         xvgr_legend(d.sfp, ngrps, grpnames, oenv);
 325     }
 326     if (fnFrac)
 327     {
 328         d.cfp = xvgropen(fnFrac, "Covered fraction", "Time (ps)", "Fraction",
 329                          oenv);
 330         xvgr_selections(d.cfp, trj);
 331         xvgr_legend(d.cfp, ngrps, grpnames, oenv);
 332     }
 333     if (fnIndex)
 334     {
 335         d.ifp = ffopen(fnIndex, "w");
 336         xvgr_selections(d.ifp, trj);
 337     }
 338     if (fnNdx)
 339     {
 340         d.block = new_blocka();
 341         d.gnames = NULL;
 342     }
 343     if (fnMask)
 344     {
 345         d.mfp = ffopen(fnMask, "w");
 346         xvgr_selections(d.mfp, trj);
 347     }
 348
 349     /* Do the analysis and write out results */
 350     gmx_ana_do(trj, 0, &print_data, &d);
 351
 352     /* Close the files */
 353     if (d.sfp)
 354     {
 355         ffclose(d.sfp);
 356     }
 357     if (d.cfp)
 358     {
 359         ffclose(d.cfp);
 360     }
 361     if (d.ifp)
 362     {
 363         ffclose(d.ifp);
 364     }
 365     if (d.block)
 366     {
 367         write_index(fnNdx, d.block, d.gnames);
 368     }
 369     if (d.mfp)
 370     {
 371         ffclose(d.mfp);
 372     }
 373
 374     thanx(stderr);
 375
 376     return 0;
 377 }