src/tools/gmx_select.c

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