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