share/template/template.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,2010,2012, 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 #include <gromacs/copyrite.h>
  39 #include <gromacs/filenm.h>
  40 #include <gromacs/macros.h>
  41 #include <gromacs/pbc.h>
  42 #include <gromacs/smalloc.h>
  43 #include <gromacs/statutil.h>
  44 #include <gromacs/xvgr.h>
  45 #include <gromacs/gmx_fatal.h>
  46
  47 #include <gromacs/nbsearch.h>
  48 #include <gromacs/trajana.h>
  49
  50 /*! \brief
  51  * Template analysis data structure.
  52  */
  53 typedef struct
  54 {
  55     gmx_ana_selection_t *refsel;
  56     FILE                *fp;
  57     real                *ave;
  58     real                *n;
  59     gmx_ana_nbsearch_t  *nb;
  60 } t_analysisdata;
  61
  62 /*! \brief
  63  * Function that does the analysis for a single frame.
  64  *
  65  * It is called once for each frame.
  66  */
  67 static int
  68 analyze_frame(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  69               int nr, gmx_ana_selection_t *sel[], void *data)
  70 {
  71     t_analysisdata     *d = (t_analysisdata *)data;
  72     int                 g, i;
  73     real                frave;
  74     int                 rc;
  75
  76     /* Here, you can do whatever analysis your program requires for a frame. */
  77     if (d->fp)
  78     {
  79         fprintf(d->fp, "%10.3f", fr->time);
  80     }
  81     rc = gmx_ana_nbsearch_pos_init(d->nb, pbc, &d->refsel->p);
  82     if (rc != 0)
  83     {
  84         gmx_fatal(FARGS, "Neighborhood search initialization failed");
  85     }
  86     for (g = 0; g < nr; ++g)
  87     {
  88         frave = 0;
  89         for (i = 0; i < sel[g]->p.nr; ++i)
  90         {
  91             frave += gmx_ana_nbsearch_pos_mindist(d->nb, &sel[g]->p, i);
  92         }
  93         d->ave[g] += frave;
  94         d->n[g]   += sel[g]->p.nr;
  95         if (d->fp)
  96         {
  97             frave /= sel[g]->p.nr;
  98             fprintf(d->fp, " %.3f", frave);
  99         }
 100     }
 101     if (d->fp)
 102     {
 103         fprintf(d->fp, "\n");
 104     }
 105     /* We need to return 0 to tell that everything went OK */
 106     return 0;
 107 }
 108
 109 /*! \brief
 110  * Function that implements the analysis tool.
 111  *
 112  * Following the style of Gromacs analysis tools, this function is called
 113  * \p gmx_something.
 114  */
 115 int
 116 gmx_template(int argc, char *argv[])
 117 {
 118     const char         *desc[] = {
 119         "This is a template for writing your own analysis tools for",
 120         "Gromacs. The advantage of using Gromacs for this is that you",
 121         "have access to all information in the topology, and your",
 122         "program will be able to handle all types of coordinates and",
 123         "trajectory files supported by Gromacs. In addition,",
 124         "you get a lot of functionality for free from the trajectory",
 125         "analysis library, including support for flexible dynamic",
 126         "selections. Go ahead an try it![PAR]",
 127         "To get started with implementing your own analysis program,",
 128         "follow the instructions in the README file provided.",
 129         "This template implements a simple analysis programs that calculates",
 130         "average distances from the a reference group to one or more",
 131         "analysis groups.",
 132     };
 133
 134     /* Command-line arguments */
 135     real                cutoff = 0;
 136     gmx_bool                bArg   = FALSE;
 137     t_pargs             pa[] = {
 138         {"-cutoff", FALSE, etREAL, {&cutoff},
 139          "Cutoff for distance calculation (0 = no cutoff)"},
 140         {"-arg2",   FALSE, etBOOL, {&bArg},
 141          "Example argument 2"},
 142     };
 143     /* The second argument is for demonstration purposes only */
 144
 145     /* Output files */
 146     t_filenm            fnm[] = {
 147         {efXVG, "-o", "avedist", ffOPTWR},
 148     };
 149 #define NFILE asize(fnm)
 150
 151     gmx_ana_traj_t       *trj;
 152     output_env_t          oenv;
 153     t_analysisdata        d;
 154     int                   ngrps;
 155     gmx_ana_selection_t **sel;
 156     int                   g;
 157     int                   rc;
 158
 159     CopyRight(stderr, argv[0]);
 160     /* Here, we can use flags to specify requirements for the selections and/or
 161      * other features of the library. */
 162     gmx_ana_traj_create(&trj, ANA_REQUIRE_TOP);
 163     gmx_ana_set_nrefgrps(trj, 1);
 164     gmx_ana_set_nanagrps(trj, -1);
 165     /* If required, other functions can also be used to configure the library
 166      * before calling parse_trjana_args(). */
 167     parse_trjana_args(trj, &argc, argv, PCA_CAN_VIEW,
 168                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
 169                       &oenv);
 170
 171     /* You can now do any initialization you wish, using the information
 172      * from the trj structure.
 173      * In particular, you should store any command-line parameter values that
 174      * the analysis part requires into d for them to be accessible. */
 175
 176     /* First, we get some selection information from the structure */
 177     gmx_ana_get_refsel(trj, 0, &d.refsel);
 178     gmx_ana_get_nanagrps(trj, &ngrps);
 179     gmx_ana_get_anagrps(trj, &sel);
 180
 181     /* First, we initialize the neighborhood search for the first index
 182      * group. */
 183     rc = gmx_ana_nbsearch_create(&d.nb, cutoff, d.refsel->p.nr);
 184     if (rc != 0)
 185     {
 186         gmx_fatal(FARGS, "neighborhood search initialization failed");
 187     }
 188
 189     /* We then allocate memory in d to store intermediate data
 190      * and initialize the counters to zero */
 191     snew(d.ave, ngrps);
 192     snew(d.n,   ngrps);
 193
 194     /* We also open the output file if the user provided it */
 195     d.fp = NULL;
 196     if (opt2bSet("-o", NFILE, fnm))
 197     {
 198         d.fp = xvgropen(opt2fn("-o", NFILE, fnm), "Average distance",
 199                         "Time [ps]", "Distance [nm]", oenv);
 200         xvgr_selections(d.fp, trj);
 201     }
 202
 203     /* Now, we do the actual analysis */
 204     gmx_ana_do(trj, 0, &analyze_frame, &d);
 205
 206     /* Now, the analysis has been done for all frames, and you can access the
 207      * results in d. Here, you should post-process your data and write out any
 208      * averaged properties. */
 209
 210     /* For the template, we close the output file if one was opened */
 211     if (d.fp)
 212     {
 213         ffclose(d.fp);
 214     }
 215
 216     /* We also calculate the average distance for each group */
 217     for (g = 0; g < ngrps; ++g)
 218     {
 219         d.ave[g] /= d.n[g];
 220         fprintf(stderr, "Average distance for '%s': %.3f nm\n",
 221                 sel[g]->name, d.ave[g]);
 222     }
 223
 224     /* Here, we could free some memory, but this usually not necessary as we
 225      * are going to quit anyways. */
 226
 227     thanx(stderr);
 228     return 0;
 229 }
 230
 231 /*! \brief
 232  * The main function.
 233  *
 234  * In Gromacs, most analysis programs are implemented such that the \p main
 235  * function is only a wrapper for a \p gmx_something function that does all
 236  * the work, and that convention is also followed here. */
 237 int
 238 main(int argc, char *argv[])
 239 {
 240     gmx_template(argc, argv);
 241     return 0;
 242 }