src/tools/g_dielectric.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_dielectric_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <stdlib.h>
  33 #include <ctype.h>
  34 #include <math.h>
  35 #include "copyrite.h"
  36 #include "typedefs.h"
  37 #include "string2.h"
  38 #include "gstat.h"
  39 #include "smalloc.h"
  40 #include "futil.h"
  41 #include "macros.h"
  42 #include "xvgr.h"
  43 #include "complex.h"
  44 #include "nr.h"
  45
  46 /* Determines at which point in the array the fit should start */
  47 int calc_nbegin(int nx,real x[],real tbegin)
  48 {
  49   int nbegin;
  50
  51   for(nbegin=0; (nbegin<nx); nbegin++) {
  52     if (x[nbegin] == tbegin)
  53       break;
  54   }
  55   fprintf(stderr,"nbegin = %d\n",nbegin);
  56
  57   return nbegin;
  58 }
  59
  60 real numerical_deriv(int nx,real x[],real y[],real fity[],real dy[],
  61                      real tendInt,int nsmooth)
  62 {
  63   FILE *tmpfp;
  64   int  i,nbegin,i0,i1;
  65   real fac,fx,fy,integralSmth;
  66   real *tmp;
  67
  68   nbegin = calc_nbegin(nx,x,tendInt);
  69   snew(tmp,nx);
  70   if (nsmooth == 0) {
  71     for(i=0; (i<nbegin); i++)
  72       tmp[i]=y[i];
  73     fac = y[nbegin]/fity[nbegin];
  74     fprintf(stderr,"scaling fitted curve by %g\n",fac);
  75     for(i=nbegin; (i<nx); i++)
  76       tmp[i]=fity[i]*fac;
  77   }
  78   else {
  79     i0 = max(0,nbegin);
  80     i1 = min(nx-1,nbegin+nsmooth);
  81     fprintf(stderr,"Making smooth transition from %d thru %d\n",i0,i1);
  82     for(i=0; (i<i0); i++)
  83       tmp[i]=y[i];
  84     for(i=i0; (i<=i1); i++) {
  85       fx = (i1-i)/(real)(i1-i0);
  86       fy = (i-i0)/(real)(i1-i0);
  87       /* fprintf(stderr,"factors for smoothing: %10g %10g\n",fx,fy);*/
  88       tmp[i] = fx*y[i] + fy*fity[i];
  89     }
  90     for(i=i1+1; (i<nx); i++)
  91       tmp[i]=fity[i];
  92   }
  93
  94   tmpfp = ffopen("integral_smth.xvg","w");
  95   integralSmth=print_and_integrate(tmpfp,nx,x[1]-x[0],tmp);
  96   fprintf(stderr,"SMOOTH integral = %10.5e\n",integralSmth);
  97
  98   dy[0] = (tmp[1]-tmp[0])/(x[1]-x[0]);
  99   for(i=1; (i<nx-1); i++) {
 100     dy[i] = (tmp[i+1]-tmp[i-1])/(x[i+1]-x[i-1]);
 101   }
 102   dy[nx-1] = (tmp[nx-1]-tmp[nx-2])/(x[nx-1]-x[nx-2]);
 103
 104   for(i=0; (i<nx); i++)
 105     dy[i] *= -1;
 106
 107   sfree(tmp);
 108
 109   return integralSmth;
 110 }
 111
 112 void do_four(char *fn,char *cn,int nx,real x[],real dy[],real eps0,real epsRF)
 113 {
 114   FILE      *fp,*cp;
 115   t_complex *tmp,gw,hw,kw;
 116   int       i,nnx,nxsav;
 117   real      fac,nu,dt,*ptr,maxeps,numax;
 118
 119   nxsav = nx;
 120   /*while ((dy[nx-1] == 0.0) && (nx > 0))
 121     nx--;*/
 122   if (nx == 0) {
 123     fprintf(stderr,"Empty dataset!\n");
 124     return;
 125   }
 126   nnx=1;
 127   while (nnx < 2*nx) {
 128     nnx*=2;
 129   }
 130   snew(tmp,2*nnx);
 131   fprintf(stderr,"Doing FFT of %d points\n",nnx);
 132   for(i=0; (i<nx); i++)
 133     tmp[i].re = dy[i];
 134   ptr=&tmp[0].re;
 135   four1(ptr-1,nnx,-1);
 136
 137   dt=x[1]-x[0];
 138   fac=((eps0-1)/(2*epsRF+eps0))/tmp[0].re;
 139   fp=xvgropen(fn,"Epsilon(\\8w\\4)","Freq. (GHz)","eps");
 140   cp=xvgropen(cn,"Cole-Cole plot","Eps'","Eps''");
 141   maxeps = 0;
 142   numax  = 0;
 143   for(i=0; (i<nxsav); i++) {
 144     gw     = rcmul(fac,tmp[i]);
 145     hw     = rcmul(2*epsRF,gw);
 146     hw.re += 1.0;
 147     gw.re  = 1.0 - gw.re;
 148     gw.im  = -gw.im;
 149     kw     = cdiv(hw,gw);
 150     kw.im *= -1;
 151
 152     nu     = (i+1)*1000.0/(nnx*dt);
 153     if (kw.im > maxeps) {
 154       maxeps = kw.im;
 155       numax  = nu;
 156     }
 157
 158     fprintf(fp,"%10.5e  %10.5e  %10.5e\n",nu,kw.re,kw.im);
 159     fprintf(cp,"%10.5e  %10.5e\n",kw.re,kw.im);
 160   }
 161   fprintf(stderr,"MAXEPS = %10.5e at frequency %10.5e GHz (tauD = %8.1f ps)\n",
 162           maxeps,numax,1000/(2*M_PI*numax));
 163   fclose(fp);
 164   fclose(cp);
 165   sfree(tmp);
 166 }
 167
 168 int main(int argc,char *argv[])
 169 {
 170   static char *desc[] = {
 171     "dielectric calculates frequency dependent dielectric constants",
 172     "from the autocorrelation function of the total dipole moment in",
 173     "your simulation. This ACF can be generated by g_dipoles.",
 174     "For an estimate of the error you can run g_statistics on the",
 175     "ACF, and use the output thus generated for this program.",
 176     "The functional forms of the available functions are:[PAR]",
 177     "One parmeter  : y = Exp[-a1 x]",
 178     "Two parmeters : y = a2 Exp[-a1 x]",
 179     "Three parmeter: y = a2 Exp[-a1 x] + (1 - a2) Exp[-a3 x]",
 180     "Startvalues for the fit procedure can be given on the commandline.",
 181     "It is also possible to fix parameters at their start value, use -nfix",
 182     "with the number of the parameter you want to fix.",
 183     "[PAR]",
 184     "Three output files are generated, the first contains the ACF,",
 185     "an exponential fit to it with 1, 2 or 3 parameters, and the",
 186     "numerical derivative of the combination data/fit.",
 187     "The second file contains the real and imaginary parts of the",
 188     "frequency-dependent dielectric constant, the last gives a plot",
 189     "known as the Cole-Cole plot, in which the  imaginary",
 190     "component is plotted as a fcuntion of the real component.",
 191     "For a pure exponential relaxation (Debye relaxation) the latter",
 192     "plot should be one half of a circle"
 193   };
 194   t_filenm fnm[] = {
 195     { efXVG, "-f", "Mtot",  ffREAD  },
 196     { efXVG, "-d", "deriv",  ffWRITE },
 197     { efXVG, "-o", "epsw",   ffWRITE },
 198     { efXVG, "-c", "cole",   ffWRITE }
 199   };
 200 #define NFILE asize(fnm)
 201   int  i,j,nx,ny,nxtail;
 202   real **y,dt,integral,fitintegral,*fitparms,fac;
 203
 204   static char *fix=NULL;
 205   static int bFour = 0,bX = 1,nfitparm=2,nsmooth=3;
 206   static real tendInt=5.0,tbegin=5.0,tend=500.0;
 207   static real A=0.0,tau1=0.0,tau2=0.0,eps0=80,epsRF=78.5,tail=500.0;
 208   real   lambda;
 209   t_pargs pa[] = {
 210     { "-fft", FALSE, etBOOL, {&bFour},
 211       "use fast fourier transform for correlation function" },
 212     { "-x1",  FALSE, etBOOL, {&bX},
 213       "use first column as X axis rather than first data set" },
 214     { "-eint", FALSE, etREAL, {&tendInt},
 215       "Time were to end the integration of the data and start to use the fit"},
 216     { "-bfit", FALSE, etREAL, {&tbegin},
 217       "Begin time of fit" },
 218     { "-efit", FALSE, etREAL, {&tend},
 219       "End time of fit" },
 220     { "-tail", FALSE, etREAL, {&tail},
 221       "Length of function including data and tail from fit" },
 222     { "-A", FALSE, etREAL, {&A},
 223       "Start value for fit parameter A" },
 224     { "-tau1", FALSE, etREAL, {&tau1},
 225       "Start value for fit parameter tau1" },
 226     { "-tau2", FALSE, etREAL, {&tau2},
 227       "Start value for fit parameter tau2" },
 228     { "-eps0", FALSE, etREAL, {&eps0},
 229       "Epsilon 0 of your liquid" },
 230     { "-epsRF", FALSE, etREAL, {&epsRF},
 231       "Epsilon of the reaction field used in your simulation" },
 232     { "-fix", FALSE, etSTR,  {&fix},
 233       "Fix this parameter at its start value, e.g. A, tau1 or tau2" },
 234     { "-nparm", FALSE, etINT,  {&nfitparm},
 235       "Number of parameters for fitting!" },
 236     { "-nsmooth", FALSE, etINT, {&nsmooth},
 237       "Number of points for smoothing" }
 238   };
 239
 240   CopyRight(stderr,argv[0]);
 241   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW,TRUE,
 242                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 243
 244   nx = read_xvg(opt2fn("-f",NFILE,fnm),&y,&ny);
 245   dt = y[0][1]-y[0][0];
 246   nxtail = max(tail/dt,nx);
 247
 248   fprintf(stderr,"Read data set containing %d colums and %d rows\n",ny,nx);
 249   fprintf(stderr,"Assuming (from data) that timestep is %g, nxtail = %d\n",
 250           dt,nxtail);
 251   if (nxtail > nx) {
 252     for(i=0; (i<ny); i++)
 253       srenew(y[i],nxtail);
 254     for(i=nx; (i<nxtail); i++) {
 255       y[0][i] = dt*i+y[0][0];
 256       for(j=1; (j<ny); j++)
 257         y[j][i] = 0.0;
 258     }
 259     nx=nxtail;
 260   }
 261
 262   /* We have read a file WITHOUT standard deviations, so we make our own... */
 263   if (ny==2) {
 264     fprintf(stderr,"Creating standard deviation numbers ...\n");
 265     srenew(y,3);
 266     snew(y[2],nx);
 267
 268     fac=1.0/((real)nx);
 269     for(i=0; (i<nx); i++)
 270       y[2][i] = fac;
 271   }
 272
 273   snew(fitparms,4);
 274   fitparms[0]=tau1;
 275   if (nfitparm > 1)
 276     fitparms[1]=A;
 277   if (nfitparm > 2)
 278     fitparms[2]=tau2;
 279
 280   if (ny < 5) {
 281     srenew(y,5);
 282     snew(y[3],nx);
 283     snew(y[4],nx);
 284   }
 285   integral = print_and_integrate(NULL,calc_nbegin(nx,y[0],tbegin),
 286                                  y[0][1]-y[0][0],y[1]);
 287   integral += do_lmfit(nx,y[1],y[2],y[0][1]-y[0][0],tbegin,tend,
 288                         TRUE,nfitparm,y[3],fitparms,fix);
 289   lambda = (eps0 - 1.0)/(2*epsRF - 1.0);
 290   fprintf(stderr,"DATA INTEGRAL: %5.1f, tauD(old) = %5.1f ps, tau_slope = %5.1f, tau_slope,D = %5.1f ps\n",
 291           integral,integral*(2*epsRF+eps0)/(2*epsRF+1),fitparms[0],
 292           fitparms[0]*(2*epsRF+eps0)/(2*epsRF+1));
 293
 294   fprintf(stderr,"tau_D from tau1 = %8.3g , eps(Infty) = %8.3f\n",
 295           fitparms[0]*(1 + fitparms[1]*lambda),
 296           1 + ((1 - fitparms[1])*(eps0 - 1))/(1 + fitparms[1]*lambda));
 297
 298   fitintegral=numerical_deriv(nx,y[0],y[1],y[3],y[4],tendInt,nsmooth);
 299   fprintf(stderr,"FIT INTEGRAL (tau_M): %5.1f, tau_D = %5.1f\n",
 300           fitintegral,fitintegral*(2*epsRF+eps0)/(2*epsRF+1));
 301   /* Now we have the negative gradient of <Phi(0) Phi(t)> */
 302
 303   dump_xvg(opt2fn("-d",NFILE,fnm),"Corr, Std Dev, Fit & Deriv",nx,5,y);
 304
 305   /* Do FFT and analysis */
 306   do_four(opt2fn("-o",NFILE,fnm),opt2fn("-c",NFILE,fnm),
 307           nx,y[0],y[4],eps0,epsRF);
 308
 309   xvgr_file(opt2fn("-o",NFILE,fnm),"-nxy -log x");
 310   xvgr_file(opt2fn("-c",NFILE,fnm),NULL);
 311   xvgr_file(opt2fn("-d",NFILE,fnm),"-nxy");
 312
 313   thanx(stdout);
 314
 315   return 0;
 316 }
 317
 318
 319
 320
 321
 322
 323
 324
 325