src/tools/acf.test.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 "typedefs.h"
  41 #include "xvgr.h"
  42 #include "gstat.h"
  43 #include "copyrite.h"
  44 #include "macros.h"
  45 #include "random.h"
  46 #include "smalloc.h"
  47
  48 int main(int argc,char *argv[])
  49 {
  50   FILE *fp;
  51   const char *desc[] = {
  52     "testac tests the functioning of the GROMACS acf routines"
  53   };
  54   static int nframes = 1024;
  55   static int datatp  = 0;
  56   static real a=0.02*M_PI;
  57   output_env_t oenv;
  58   t_pargs pa[] = {
  59     { "-np", FALSE, etINT, &nframes,
  60       "Number of data points" },
  61     { "-dtp",FALSE, etINT, &datatp,
  62       "Which data: 0=all 0.0, 1=all 1.0, 2=cos(a t), 3=random, 4=cos(a t)+random, 5=sin(a t)/(a t)" }
  63   };
  64   static char *str[] = {
  65     "all 0.0",
  66     "all 1.0",
  67     "cos(a t)",
  68     "random",
  69     "cos(a t)+random",
  70     "sin(a t)/(a t)"
  71   };
  72   t_filenm fnm[] = {
  73     { efXVG, "-d", "acf-data", ffWRITE },
  74     { efXVG, "-c", "acf-corr", ffWRITE },
  75     { efXVG, "-comb", "acf-comb.xvg", ffWRITE }
  76   };
  77 #define NFILE asize(fnm)
  78   int     npargs,i,nlag;
  79   int     seed=1993;
  80   real    *data,*data2,x;
  81   t_pargs *ppa;
  82
  83   CopyRight(stderr,argv[0]);
  84   npargs = asize(pa);
  85   ppa    = add_acf_pargs(&npargs,pa);
  86   parse_common_args_r(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
  87                       NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
  88   snew(data,nframes);
  89   snew(data2,nframes);
  90
  91   fp = xvgropen(opt2fn("-d",NFILE,fnm),"testac","x","y",oenv);
  92   for(i=0; (i<nframes); i++) {
  93     x = a*i;
  94     switch (datatp) {
  95     case 1:
  96       data[i] = 1;
  97       break;
  98     case 2:
  99       data[i] = cos(x);
 100       break;
 101     case 3:
 102       data[i] = 2*rando(&seed)-1.0;
 103       break;
 104     case 4:
 105       data[i] = cos(x)+2*rando(&seed)-1.0;
 106       break;
 107     case 5:
 108       if (i==0)
 109         data[i] = 1;
 110       else
 111         data[i] = sin(x)/(x);
 112     default:
 113       /* Data remains 0.0 */
 114       break;
 115     }
 116     fprintf(fp,"%10g  %10g\n",x,data[i]);
 117     data2[i] = data[i];
 118   }
 119   fclose(fp);
 120
 121   do_autocorr(opt2fn("-c",NFILE,fnm),oenv,str[datatp],
 122               nframes,1,&data,a,eacNormal,FALSE);
 123
 124   nlag = get_acfnout();
 125   fp = xvgropen(opt2fn("-comb",NFILE,fnm),"testac","x","y",oenv);
 126   for(i=0; (i<nlag); i++) {
 127     fprintf(fp,"%10g  %10g  %10g\n",a*i,data2[i],data[i]);
 128   }
 129   ffclose(fp);
 130
 131   do_view(opt2fn("-c",NFILE,fnm),"-nxy");
 132
 133   thanx(stderr);
 134
 135   return 0;
 136 }