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     gmx_bool                bDump;
  56     gmx_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         "[TT]g_select[tt] 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.",
 184         "[TT]-os[tt] is the default output option if none is selected.[PAR]",
 185         "With [TT]-os[tt], calculates the number of positions in each",
 186         "selection for each frame. With [TT]-norm[tt], the output is",
 187         "between 0 and 1 and describes the fraction from the maximum",
 188         "number of positions (e.g., for selection 'resname RA and x < 5'",
 189         "the maximum number of positions is the number of atoms in",
 190         "RA residues). With [TT]-cfnorm[tt], the output is divided",
 191         "by the fraction covered by the selection.",
 192         "[TT]-norm[tt] and [TT]-cfnorm[tt] can be specified independently",
 193         "of one another.[PAR]",
 194         "With [TT]-oc[tt], the fraction covered by each selection is",
 195         "written out as a function of time.[PAR]",
 196         "With [TT]-oi[tt], the selected atoms/residues/molecules are",
 197         "written out as a function of time. In the output, the first",
 198         "column contains the frame time, the second contains the number",
 199         "of positions, followed by the atom/residue/molecule numbers.",
 200         "If more than one selection is specified, the size of the second",
 201         "group immediately follows the last number of the first group",
 202         "and so on. With [TT]-dump[tt], the frame time and the number",
 203         "of positions is omitted from the output. In this case, only one",
 204         "selection can be given.[PAR]",
 205         "With [TT]-on[tt], the selected atoms are written as a index file",
 206         "compatible with [TT]make_ndx[tt] and the analyzing tools. Each selection",
 207         "is written as a selection group and for dynamic selections a",
 208         "group is written for each frame.[PAR]",
 209         "For residue numbers, the output of [TT]-oi[tt] can be controlled",
 210         "with [TT]-resnr[tt]: [TT]number[tt] (default) prints the residue",
 211         "numbers as they appear in the input file, while [TT]index[tt] prints",
 212         "unique numbers assigned to the residues in the order they appear",
 213         "in the input file, starting with 1. The former is more intuitive,",
 214         "but if the input contains multiple residues with the same number,",
 215         "the output can be less useful.[PAR]",
 216         "With [TT]-om[tt], a mask is printed for the first selection",
 217         "as a function of time. Each line in the output corresponds to",
 218         "one frame, and contains either 0/1 for each atom/residue/molecule",
 219         "possibly selected. 1 stands for the atom/residue/molecule being",
 220         "selected for the current frame, 0 for not selected.",
 221         "With [TT]-dump[tt], the frame time is omitted from the output.",
 222     };
 223
 224     gmx_bool                bDump     = FALSE;
 225     gmx_bool                bFracNorm = FALSE;
 226     gmx_bool                bTotNorm  = FALSE;
 227     const char         *routt[] = {NULL, "number", "index", NULL};
 228     t_pargs             pa[] = {
 229         {"-dump",   FALSE, etBOOL, {&bDump},
 230          "Do not print the frame time (-om, -oi) or the index size (-oi)"},
 231         {"-norm",   FALSE, etBOOL, {&bTotNorm},
 232          "Normalize by total number of positions with -os"},
 233         {"-cfnorm", FALSE, etBOOL, {&bFracNorm},
 234          "Normalize by covered fraction with -os"},
 235         {"-resnr",  FALSE, etENUM, {routt},
 236          "Residue number output type"},
 237     };
 238
 239     t_filenm            fnm[] = {
 240         {efXVG, "-os", "size.xvg",  ffOPTWR},
 241         {efXVG, "-oc", "cfrac.xvg", ffOPTWR},
 242         {efDAT, "-oi", "index.dat", ffOPTWR},
 243         {efDAT, "-om", "mask.dat",  ffOPTWR},
 244         {efNDX, "-on", "index.ndx", ffOPTWR},
 245     };
 246 #define NFILE asize(fnm)
 247
 248     gmx_ana_traj_t       *trj;
 249     t_topology           *top;
 250     int                   ngrps;
 251     gmx_ana_selection_t **sel;
 252     char                **grpnames;
 253     t_dsdata              d;
 254     const char            *fnSize, *fnFrac, *fnIndex, *fnNdx, *fnMask;
 255     int                   g;
 256     int                   rc;
 257     output_env_t          oenv;
 258
 259     CopyRight(stderr, argv[0]);
 260     gmx_ana_traj_create(&trj, 0);
 261     gmx_ana_set_nanagrps(trj, -1);
 262     parse_trjana_args(trj, &argc, argv, 0,
 263                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
 264                       &oenv);
 265     gmx_ana_get_nanagrps(trj, &ngrps);
 266     gmx_ana_get_anagrps(trj, &sel);
 267     gmx_ana_init_coverfrac(trj, CFRAC_SOLIDANGLE);
 268
 269     /* Get output file names */
 270     fnSize  = opt2fn_null("-os", NFILE, fnm);
 271     fnFrac  = opt2fn_null("-oc", NFILE, fnm);
 272     fnIndex = opt2fn_null("-oi", NFILE, fnm);
 273     fnNdx   = opt2fn_null("-on", NFILE, fnm);
 274     fnMask  = opt2fn_null("-om", NFILE, fnm);
 275     /* Write out sizes if nothing specified */
 276     if (!fnFrac && !fnIndex && !fnMask && !fnNdx)
 277     {
 278         fnSize = opt2fn("-os", NFILE, fnm);
 279     }
 280
 281     if ( bDump && ngrps > 1)
 282     {
 283         gmx_fatal(FARGS, "Only one index group allowed with -dump");
 284     }
 285     if (fnNdx && sel[0]->p.m.type != INDEX_ATOM)
 286     {
 287         gmx_fatal(FARGS, "Only atom selection allowed with -on");
 288     }
 289     if (fnMask && ngrps > 1)
 290     {
 291         fprintf(stderr, "warning: the mask (-om) will only be written for the first group\n");
 292     }
 293     if (fnMask && !sel[0]->bDynamic)
 294     {
 295         fprintf(stderr, "warning: will not write the mask (-om) for a static selection\n");
 296         fnMask = NULL;
 297     }
 298
 299     /* Initialize reference calculation for masks */
 300     if (fnMask)
 301     {
 302         gmx_ana_get_topology(trj, FALSE, &top, NULL);
 303         snew(d.mmap, 1);
 304         gmx_ana_indexmap_init(d.mmap, sel[0]->g, top, sel[0]->p.m.type);
 305     }
 306
 307     /* Initialize calculation data */
 308     d.bDump     = bDump;
 309     d.bFracNorm = bFracNorm;
 310     d.routt     = routt[0];
 311     snew(d.size,  ngrps);
 312     for (g = 0; g < ngrps; ++g)
 313     {
 314         d.size[g] = bTotNorm ? sel[g]->p.nr : 1;
 315     }
 316
 317     /* Open output files */
 318     d.sfp = d.cfp = d.ifp = d.mfp = NULL;
 319     d.block = NULL;
 320     gmx_ana_get_grpnames(trj, &grpnames);
 321     if (fnSize)
 322     {
 323         d.sfp = xvgropen(fnSize, "Selection size", "Time (ps)", "Number",oenv);
 324         xvgr_selections(d.sfp, trj);
 325         xvgr_legend(d.sfp, ngrps, (const char**)grpnames, oenv);
 326     }
 327     if (fnFrac)
 328     {
 329         d.cfp = xvgropen(fnFrac, "Covered fraction", "Time (ps)", "Fraction",
 330                          oenv);
 331         xvgr_selections(d.cfp, trj);
 332         xvgr_legend(d.cfp, ngrps, (const char**)grpnames, oenv);
 333     }
 334     if (fnIndex)
 335     {
 336         d.ifp = ffopen(fnIndex, "w");
 337         xvgr_selections(d.ifp, trj);
 338     }
 339     if (fnNdx)
 340     {
 341         d.block = new_blocka();
 342         d.gnames = NULL;
 343     }
 344     if (fnMask)
 345     {
 346         d.mfp = ffopen(fnMask, "w");
 347         xvgr_selections(d.mfp, trj);
 348     }
 349
 350     /* Do the analysis and write out results */
 351     gmx_ana_do(trj, 0, &print_data, &d);
 352
 353     /* Close the files */
 354     if (d.sfp)
 355     {
 356         ffclose(d.sfp);
 357     }
 358     if (d.cfp)
 359     {
 360         ffclose(d.cfp);
 361     }
 362     if (d.ifp)
 363     {
 364         ffclose(d.ifp);
 365     }
 366     if (d.block)
 367     {
 368         write_index(fnNdx, d.block, d.gnames);
 369     }
 370     if (d.mfp)
 371     {
 372         ffclose(d.mfp);
 373     }
 374
 375     thanx(stderr);
 376
 377     return 0;
 378 }