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