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