src/gromacs/gmxana/gmx_rotacf.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  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <string.h>
  41 #include "sysstuff.h"
  42 #include "physics.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "futil.h"
  46 #include "statutil.h"
  47 #include "index.h"
  48 #include "macros.h"
  49 #include "gmx_fatal.h"
  50 #include "xvgr.h"
  51 #include "gstat.h"
  52 #include "vec.h"
  53 #include "gmx_ana.h"
  54
  55
  56 int gmx_rotacf(int argc, char *argv[])
  57 {
  58     const char     *desc[] = {
  59         "[TT]g_rotacf[tt] calculates the rotational correlation function",
  60         "for molecules. Atom triplets (i,j,k) must be given in the index",
  61         "file, defining two vectors ij and jk. The rotational ACF",
  62         "is calculated as the autocorrelation function of the vector",
  63         "n = ij x jk, i.e. the cross product of the two vectors.",
  64         "Since three atoms span a plane, the order of the three atoms",
  65         "does not matter. Optionally, by invoking the [TT]-d[tt] switch, you can",
  66         "calculate the rotational correlation function for linear molecules",
  67         "by specifying atom pairs (i,j) in the index file.",
  68         "[PAR]",
  69         "EXAMPLES[PAR]",
  70         "[TT]g_rotacf -P 1 -nparm 2 -fft -n index -o rotacf-x-P1",
  71         "-fa expfit-x-P1 -beginfit 2.5 -endfit 20.0[tt][PAR]",
  72         "This will calculate the rotational correlation function using a first",
  73         "order Legendre polynomial of the angle of a vector defined by the index",
  74         "file. The correlation function will be fitted from 2.5 ps until 20.0 ps",
  75         "to a two-parameter exponential."
  76     };
  77     static gmx_bool bVec    = FALSE, bAver = TRUE;
  78
  79     t_pargs         pa[] = {
  80         { "-d",   FALSE, etBOOL, {&bVec},
  81           "Use index doublets (vectors) for correlation function instead of triplets (planes)" },
  82         { "-aver", FALSE, etBOOL, {&bAver},
  83           "Average over molecules" }
  84     };
  85
  86     t_trxstatus    *status;
  87     int             isize;
  88     atom_id        *index;
  89     char           *grpname;
  90     rvec           *x, *x_s;
  91     matrix          box;
  92     real          **c1;
  93     rvec            xij, xjk, n;
  94     int             i, m, teller, n_alloc, natoms, nvec, ai, aj, ak;
  95     unsigned long   mode;
  96     real            t, t0, t1, dt;
  97     gmx_rmpbc_t     gpbc = NULL;
  98     t_topology     *top;
  99     int             ePBC;
 100     t_filenm        fnm[] = {
 101         { efTRX, "-f", NULL,  ffREAD  },
 102         { efTPX, NULL, NULL,  ffREAD },
 103         { efNDX, NULL, NULL,  ffREAD  },
 104         { efXVG, "-o", "rotacf",  ffWRITE }
 105     };
 106 #define NFILE asize(fnm)
 107     int             npargs;
 108     t_pargs        *ppa;
 109
 110     output_env_t    oenv;
 111
 112     npargs = asize(pa);
 113     ppa    = add_acf_pargs(&npargs, pa);
 114
 115     parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 116                       NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv);
 117
 118     rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &isize, &index, &grpname);
 119
 120     if (bVec)
 121     {
 122         nvec = isize/2;
 123     }
 124     else
 125     {
 126         nvec = isize/3;
 127     }
 128
 129     if (((isize % 3) != 0) && !bVec)
 130     {
 131         gmx_fatal(FARGS, "number of index elements not multiple of 3, "
 132                   "these can not be atom triplets\n");
 133     }
 134     if (((isize % 2) != 0) && bVec)
 135     {
 136         gmx_fatal(FARGS, "number of index elements not multiple of 2, "
 137                   "these can not be atom doublets\n");
 138     }
 139
 140     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
 141
 142     snew(c1, nvec);
 143     for (i = 0; (i < nvec); i++)
 144     {
 145         c1[i] = NULL;
 146     }
 147     n_alloc = 0;
 148
 149     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 150     snew(x_s, natoms);
 151
 152     gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms);
 153
 154     /* Start the loop over frames */
 155     t1      = t0 = t;
 156     teller  = 0;
 157     do
 158     {
 159         if (teller >= n_alloc)
 160         {
 161             n_alloc += 100;
 162             for (i = 0; (i < nvec); i++)
 163             {
 164                 srenew(c1[i], DIM*n_alloc);
 165             }
 166         }
 167         t1 = t;
 168
 169         /* Remove periodicity */
 170         gmx_rmpbc_copy(gpbc, natoms, box, x, x_s);
 171
 172         /* Compute crossproducts for all vectors, if triplets.
 173          * else, just get the vectors in case of doublets.
 174          */
 175         if (bVec == FALSE)
 176         {
 177             for (i = 0; (i < nvec); i++)
 178             {
 179                 ai = index[3*i];
 180                 aj = index[3*i+1];
 181                 ak = index[3*i+2];
 182                 rvec_sub(x_s[ai], x_s[aj], xij);
 183                 rvec_sub(x_s[aj], x_s[ak], xjk);
 184                 cprod(xij, xjk, n);
 185                 for (m = 0; (m < DIM); m++)
 186                 {
 187                     c1[i][DIM*teller+m] = n[m];
 188                 }
 189             }
 190         }
 191         else
 192         {
 193             for (i = 0; (i < nvec); i++)
 194             {
 195                 ai = index[2*i];
 196                 aj = index[2*i+1];
 197                 rvec_sub(x_s[ai], x_s[aj], n);
 198                 for (m = 0; (m < DIM); m++)
 199                 {
 200                     c1[i][DIM*teller+m] = n[m];
 201                 }
 202             }
 203         }
 204         /* Increment loop counter */
 205         teller++;
 206     }
 207     while (read_next_x(oenv, status, &t, x, box));
 208     close_trj(status);
 209     fprintf(stderr, "\nDone with trajectory\n");
 210
 211     gmx_rmpbc_done(gpbc);
 212
 213
 214     /* Autocorrelation function */
 215     if (teller < 2)
 216     {
 217         fprintf(stderr, "Not enough frames for correlation function\n");
 218     }
 219     else
 220     {
 221         dt = (t1 - t0)/(teller-1);
 222
 223         mode = eacVector;
 224
 225         do_autocorr(ftp2fn(efXVG, NFILE, fnm), oenv, "Rotational Correlation Function",
 226                     teller, nvec, c1, dt, mode, bAver);
 227     }
 228
 229     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), NULL);
 230
 231     return 0;
 232 }