src/tools/g_rotacf.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_rotacf_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "sysstuff.h"
  34 #include "physics.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "futil.h"
  38 #include "statutil.h"
  39 #include "copyrite.h"
  40 #include "rdgroup.h"
  41 #include "macros.h"
  42 #include "fatal.h"
  43 #include "xvgr.h"
  44 #include "gstat.h"
  45 #include "vec.h"
  46
  47 int main(int argc,char *argv[])
  48 {
  49   static char *desc[] = {
  50     "g_rotacf calculates the rotational correlation function",
  51     "for molecules. Three atoms (i,j,k) must be given in the index",
  52     "file, defining two vectors ij and jk. The rotational acf",
  53     "is calculated as the autocorrelation function of the vector",
  54     "n = ij x jk, i.e. the cross product of the two vectors.",
  55     "Since three atoms span a plane, the order of the three atoms",
  56     "does not matter. Optionally, controlled by the -d switch, you can",
  57     "calculate the rotational correlation function for linear molecules",
  58     "by specifying two atoms (i,j) in the index file.",
  59     "[PAR]",
  60     "EXAMPLES[PAR]",
  61     "g_rotacf -P 1 -nparm 2 -fft -n index -o rotacf-x-P1",
  62     "-fa expfit-x-P1 -beginfit 2.5 -endfit 20.0[PAR]",
  63     "This will calculate the rotational correlation function using a first",
  64     "order Legendre polynomial of the angle of a vector defined by the index",
  65     "file. The correlation function will be fitted from 2.5 ps till 20.0 ps",
  66     "to a two parameter exponential",
  67
  68
  69     ""
  70   };
  71   static bool bVec    = FALSE;
  72
  73   t_pargs pa[] = {
  74     { "-d",   FALSE, etBOOL, {&bVec},
  75       "Use index doublets (vectors) for correlation function instead of triplets (planes)" }
  76   };
  77
  78   int        status,isize;
  79   atom_id    *index;
  80   char       *grpname;
  81   rvec       *x,*x_s;
  82   matrix     box;
  83   real       **c1;
  84   rvec       xij,xjk,n;
  85   int        i,m,teller,n_alloc,natoms,nvec,ai,aj,ak;
  86   unsigned long mode;
  87   real       t,t0,t1,dt;
  88   t_topology *top;
  89   t_filenm   fnm[] = {
  90     { efTRX, "-f", NULL,  ffREAD  },
  91     { efTPX, NULL, NULL,  ffREAD },
  92     { efNDX, NULL, NULL,  ffREAD  },
  93     { efXVG, "-o", "rotacf",  ffWRITE }
  94   };
  95 #define NFILE asize(fnm)
  96   int     npargs;
  97   t_pargs *ppa;
  98
  99   CopyRight(stderr,argv[0]);
 100   npargs = asize(pa);
 101   ppa    = add_acf_pargs(&npargs,pa);
 102
 103   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 104                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL);
 105
 106   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize,&index,&grpname);
 107
 108   if (bVec)
 109     nvec = isize/2;
 110   else
 111     nvec = isize/3;
 112
 113   if (((isize % 3) != 0) && !bVec)
 114     fatal_error(0,"number of index elements not multiple of 3, "
 115                 "these can not be atom triplets\n");
 116   if (((isize % 2) != 0) && bVec)
 117     fatal_error(0,"number of index elements not multiple of 2, "
 118                 "these can not be atom doublets\n");
 119
 120   top=read_top(ftp2fn(efTPX,NFILE,fnm));
 121
 122   snew(c1,nvec);
 123   for (i=0; (i<nvec); i++)
 124     c1[i]=NULL;
 125   n_alloc=0;
 126
 127   natoms=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 128   snew(x_s,natoms);
 129
 130   /* Start the loop over frames */
 131   t1 = t0 = t;
 132   teller  = 0;
 133   do {
 134     if (teller >= n_alloc) {
 135       n_alloc+=100;
 136       for (i=0; (i<nvec); i++)
 137         srenew(c1[i],DIM*n_alloc);
 138     }
 139     t1 = t;
 140
 141     /* Remove periodicity */
 142     rm_pbc(&(top->idef),natoms,box,x,x_s);
 143
 144     /* Compute crossproducts for all vectors, if triplets.
 145      * else, just get the vectors in case of doublets.
 146      */
 147     if (bVec == FALSE) {
 148       for (i=0; (i<nvec); i++) {
 149         ai=index[3*i];
 150         aj=index[3*i+1];
 151         ak=index[3*i+2];
 152         rvec_sub(x_s[ai],x_s[aj],xij);
 153         rvec_sub(x_s[aj],x_s[ak],xjk);
 154         oprod(xij,xjk,n);
 155         for(m=0; (m<DIM); m++)
 156           c1[i][DIM*teller+m]=n[m];
 157       }
 158     }
 159     else {
 160       for (i=0; (i<nvec); i++) {
 161         ai=index[2*i];
 162         aj=index[2*i+1];
 163         rvec_sub(x_s[ai],x_s[aj],n);
 164         for(m=0; (m<DIM); m++)
 165           c1[i][DIM*teller+m]=n[m];
 166       }
 167     }
 168     /* Increment loop counter */
 169     teller++;
 170   } while (read_next_x(status,&t,natoms,x,box));
 171   close_trj(status);
 172   fprintf(stderr,"\nDone with trajectory\n");
 173
 174   /* Autocorrelation function */
 175   if (teller < 2)
 176     fprintf(stderr,"Not enough frames for correlation function\n");
 177   else {
 178     dt=(t1 - t0)/(teller-1);
 179
 180     mode = eacVector;
 181
 182     do_autocorr(ftp2fn(efXVG,NFILE,fnm),"Rotational Correlation Function",
 183                 teller,nvec,c1,dt,mode,TRUE);
 184   }
 185
 186   xvgr_file(ftp2fn(efXVG,NFILE,fnm),NULL);
 187
 188   thanx(stdout);
 189
 190   return 0;
 191 }