src/tools/gmx_dielectric.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <stdio.h>
  42 #include <stdlib.h>
  43 #include <ctype.h>
  44 #include <math.h>
  45
  46 #include "copyrite.h"
  47 #include "typedefs.h"
  48 #include "string2.h"
  49 #include "gstat.h"
  50 #include "smalloc.h"
  51 #include "futil.h"
  52 #include "macros.h"
  53 #include "maths.h"
  54 #include "xvgr.h"
  55 #include "gmxcomplex.h"
  56 #include "correl.h"
  57 #include "gmx_ana.h"
  58 #include "gmx_fatal.h"
  59
  60
  61 /* Determines at which point in the array the fit should start */
  62 int calc_nbegin(int nx,real x[],real tbegin)
  63 {
  64   int  nbegin;
  65   real dt,dtt;
  66
  67   /* Assume input x is sorted */
  68   for(nbegin=0; (nbegin < nx) && (x[nbegin] <= tbegin); nbegin++)
  69     ;
  70   if ((nbegin == nx) || (nbegin == 0))
  71     gmx_fatal(FARGS,"Begin time %f not in x-domain [%f through %f]\n",
  72                 tbegin,x[0],x[nx-1]);
  73
  74   /* Take the one closest to tbegin */
  75   if (fabs(x[nbegin]-tbegin) > fabs(x[nbegin-1]-tbegin))
  76     nbegin--;
  77
  78   printf("nbegin = %d, x[nbegin] = %g, tbegin = %g\n",
  79           nbegin,x[nbegin],tbegin);
  80
  81   return nbegin;
  82 }
  83
  84 real numerical_deriv(int nx,real x[],real y[],real fity[],real combined[],real dy[],
  85                      real tendInt,int nsmooth)
  86 {
  87   FILE *tmpfp;
  88   int  i,nbegin,i0,i1;
  89   real fac,fx,fy,integralSmth;
  90
  91   nbegin = calc_nbegin(nx,x,tendInt);
  92   if (nsmooth == 0) {
  93     for(i=0; (i<nbegin); i++)
  94       combined[i]=y[i];
  95     fac = y[nbegin]/fity[nbegin];
  96     printf("scaling fitted curve by %g\n",fac);
  97     for(i=nbegin; (i<nx); i++)
  98       combined[i]=fity[i]*fac;
  99   }
 100   else {
 101     i0 = max(0,nbegin);
 102     i1 = min(nx-1,nbegin+nsmooth);
 103     printf("Making smooth transition from %d through %d\n",i0,i1);
 104     for(i=0; (i<i0); i++)
 105       combined[i]=y[i];
 106     for(i=i0; (i<=i1); i++) {
 107       fx = (i1-i)/(real)(i1-i0);
 108       fy = (i-i0)/(real)(i1-i0);
 109       if (debug)
 110         fprintf(debug,"x: %g factors for smoothing: %10g %10g\n",x[i],fx,fy);
 111       combined[i] = fx*y[i] + fy*fity[i];
 112     }
 113     for(i=i1+1; (i<nx); i++)
 114       combined[i]=fity[i];
 115   }
 116
 117   tmpfp = ffopen("integral_smth.xvg","w");
 118   integralSmth=print_and_integrate(tmpfp,nx,x[1]-x[0],combined,NULL,1);
 119   printf("SMOOTH integral = %10.5e\n",integralSmth);
 120
 121   dy[0] = (combined[1]-combined[0])/(x[1]-x[0]);
 122   for(i=1; (i<nx-1); i++) {
 123     dy[i] = (combined[i+1]-combined[i-1])/(x[i+1]-x[i-1]);
 124   }
 125   dy[nx-1] = (combined[nx-1]-combined[nx-2])/(x[nx-1]-x[nx-2]);
 126
 127   for(i=0; (i<nx); i++)
 128     dy[i] *= -1;
 129
 130   return integralSmth;
 131 }
 132
 133 void do_four(const char *fn,const char *cn,int nx,real x[],real dy[],
 134              real eps0,real epsRF, const output_env_t oenv)
 135 {
 136   FILE      *fp,*cp;
 137   t_complex *tmp,gw,hw,kw;
 138   int       i,nnx,nxsav;
 139   real      fac,nu,dt,*ptr,maxeps,numax;
 140
 141   nxsav = nx;
 142   /*while ((dy[nx-1] == 0.0) && (nx > 0))
 143     nx--;*/
 144   if (nx == 0)
 145     gmx_fatal(FARGS,"Empty dataset in %s, line %d!",__FILE__,__LINE__);
 146
 147   nnx=1;
 148   while (nnx < 2*nx) {
 149     nnx*=2;
 150   }
 151   snew(tmp,2*nnx);
 152   printf("Doing FFT of %d points\n",nnx);
 153   for(i=0; (i<nx); i++)
 154     tmp[i].re = dy[i];
 155   ptr=&tmp[0].re;
 156   four1(ptr-1,nnx,-1);
 157
 158   dt=x[1]-x[0];
 159   if (epsRF == 0)
 160     fac = (eps0-1)/tmp[0].re;
 161   else
 162     fac=((eps0-1)/(2*epsRF+eps0))/tmp[0].re;
 163   fp=xvgropen(fn,"Epsilon(\\8w\\4)","Freq. (GHz)","eps",oenv);
 164   cp=xvgropen(cn,"Cole-Cole plot","Eps'","Eps''",oenv);
 165   maxeps = 0;
 166   numax  = 0;
 167   for(i=0; (i<nxsav); i++) {
 168     if (epsRF == 0) {
 169       kw.re = 1+fac*tmp[i].re;
 170       kw.im = 1+fac*tmp[i].im;
 171     }
 172     else {
 173       gw     = rcmul(fac,tmp[i]);
 174       hw     = rcmul(2*epsRF,gw);
 175       hw.re += 1.0;
 176       gw.re  = 1.0 - gw.re;
 177       gw.im  = -gw.im;
 178       kw     = cdiv(hw,gw);
 179     }
 180     kw.im *= -1;
 181
 182     nu     = (i+1)*1000.0/(nnx*dt);
 183     if (kw.im > maxeps) {
 184       maxeps = kw.im;
 185       numax  = nu;
 186     }
 187
 188     fprintf(fp,"%10.5e  %10.5e  %10.5e\n",nu,kw.re,kw.im);
 189     fprintf(cp,"%10.5e  %10.5e\n",kw.re,kw.im);
 190   }
 191   printf("MAXEPS = %10.5e at frequency %10.5e GHz (tauD = %8.1f ps)\n",
 192           maxeps,numax,1000/(2*M_PI*numax));
 193   ffclose(fp);
 194   ffclose(cp);
 195   sfree(tmp);
 196 }
 197
 198 int gmx_dielectric(int argc,char *argv[])
 199 {
 200   const char *desc[] = {
 201     "[TT]g_dielectric[tt] calculates frequency dependent dielectric constants",
 202     "from the autocorrelation function of the total dipole moment in",
 203     "your simulation. This ACF can be generated by [TT]g_dipoles[tt].",
 204     "The functional forms of the available functions are:[PAR]",
 205     "One parameter:    y = [EXP]-a[SUB]1[sub] x[exp],[BR]",
 206     "Two parameters:   y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp],[BR]",
 207     "Three parameters: y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp] + (1 - a[SUB]2[sub]) [EXP]-a[SUB]3[sub] x[exp].[BR]",
 208     "Start values for the fit procedure can be given on the command line.",
 209     "It is also possible to fix parameters at their start value, use [TT]-fix[tt]",
 210     "with the number of the parameter you want to fix.",
 211     "[PAR]",
 212     "Three output files are generated, the first contains the ACF,",
 213     "an exponential fit to it with 1, 2 or 3 parameters, and the",
 214     "numerical derivative of the combination data/fit.",
 215     "The second file contains the real and imaginary parts of the",
 216     "frequency-dependent dielectric constant, the last gives a plot",
 217     "known as the Cole-Cole plot, in which the imaginary",
 218     "component is plotted as a function of the real component.",
 219     "For a pure exponential relaxation (Debye relaxation) the latter",
 220     "plot should be one half of a circle."
 221   };
 222   t_filenm fnm[] = {
 223     { efXVG, "-f", "dipcorr",ffREAD  },
 224     { efXVG, "-d", "deriv",  ffWRITE },
 225     { efXVG, "-o", "epsw",   ffWRITE },
 226     { efXVG, "-c", "cole",   ffWRITE }
 227   };
 228 #define NFILE asize(fnm)
 229   output_env_t oenv;
 230   int  i,j,nx,ny,nxtail,eFitFn,nfitparm;
 231   real dt,integral,fitintegral,*fitparms,fac,rffac;
 232   double **yd;
 233   real   **y;
 234   const char *legend[] = { "Correlation", "Std. Dev.", "Fit", "Combined", "Derivative" };
 235   static int fix=0,bFour = 0,bX = 1,nsmooth=3;
 236   static real tendInt=5.0,tbegin=5.0,tend=500.0;
 237   static real A=0.5,tau1=10.0,tau2=1.0,eps0=80,epsRF=78.5,tail=500.0;
 238   real   lambda;
 239   t_pargs pa[] = {
 240     { "-fft", FALSE, etBOOL, {&bFour},
 241       "use fast fourier transform for correlation function" },
 242     { "-x1",  FALSE, etBOOL, {&bX},
 243       "use first column as [IT]x[it]-axis rather than first data set" },
 244     { "-eint", FALSE, etREAL, {&tendInt},
 245       "Time to end the integration of the data and start to use the fit"},
 246     { "-bfit", FALSE, etREAL, {&tbegin},
 247       "Begin time of fit" },
 248     { "-efit", FALSE, etREAL, {&tend},
 249       "End time of fit" },
 250     { "-tail", FALSE, etREAL, {&tail},
 251       "Length of function including data and tail from fit" },
 252     { "-A", FALSE, etREAL, {&A},
 253       "Start value for fit parameter A" },
 254     { "-tau1", FALSE, etREAL, {&tau1},
 255       "Start value for fit parameter [GRK]tau[grk]1" },
 256     { "-tau2", FALSE, etREAL, {&tau2},
 257       "Start value for fit parameter [GRK]tau[grk]2" },
 258     { "-eps0", FALSE, etREAL, {&eps0},
 259       "[GRK]epsilon[grk]0 of your liquid" },
 260     { "-epsRF", FALSE, etREAL, {&epsRF},
 261       "[GRK]epsilon[grk] of the reaction field used in your simulation. A value of 0 means infinity." },
 262     { "-fix", FALSE, etINT,  {&fix},
 263       "Fix parameters at their start values, A (2), tau1 (1), or tau2 (4)" },
 264     { "-ffn",    FALSE, etENUM, {s_ffn},
 265       "Fit function" },
 266     { "-nsmooth", FALSE, etINT, {&nsmooth},
 267       "Number of points for smoothing" }
 268   };
 269
 270   CopyRight(stderr,argv[0]);
 271   parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_CAN_VIEW | PCA_BE_NICE,
 272                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 273   please_cite(stdout,"Spoel98a");
 274   printf("WARNING: non-polarizable models can never yield an infinite\n"
 275          "dielectric constant that is different from 1. This is incorrect\n"
 276          "in the reference given above (Spoel98a).\n\n");
 277
 278
 279   nx     = read_xvg(opt2fn("-f",NFILE,fnm),&yd,&ny);
 280   dt     = yd[0][1] - yd[0][0];
 281   nxtail = min(tail/dt,nx);
 282
 283   printf("Read data set containing %d colums and %d rows\n",ny,nx);
 284   printf("Assuming (from data) that timestep is %g, nxtail = %d\n",
 285           dt,nxtail);
 286   snew(y,6);
 287   for(i=0; (i<ny); i++)
 288     snew(y[i],max(nx,nxtail));
 289   for(i=0; (i<nx); i++) {
 290     y[0][i] = yd[0][i];
 291     for(j=1; (j<ny); j++)
 292       y[j][i] = yd[j][i];
 293   }
 294   if (nxtail > nx) {
 295     for(i=nx; (i<nxtail); i++) {
 296       y[0][i] = dt*i+y[0][0];
 297       for(j=1; (j<ny); j++)
 298         y[j][i] = 0.0;
 299     }
 300     nx=nxtail;
 301   }
 302
 303
 304   /* We have read a file WITHOUT standard deviations, so we make our own... */
 305   if (ny==2) {
 306     printf("Creating standard deviation numbers ...\n");
 307     srenew(y,3);
 308     snew(y[2],nx);
 309
 310     fac=1.0/((real)nx);
 311     for(i=0; (i<nx); i++)
 312       y[2][i] = fac;
 313   }
 314
 315   eFitFn = sffn2effn(s_ffn);
 316   nfitparm = nfp_ffn[eFitFn];
 317   snew(fitparms,4);
 318   fitparms[0]=tau1;
 319   if (nfitparm > 1)
 320     fitparms[1]=A;
 321   if (nfitparm > 2)
 322     fitparms[2]=tau2;
 323
 324
 325   snew(y[3],nx);
 326   snew(y[4],nx);
 327   snew(y[5],nx);
 328
 329   integral = print_and_integrate(NULL,calc_nbegin(nx,y[0],tbegin),
 330                                  dt,y[1],NULL,1);
 331   integral += do_lmfit(nx,y[1],y[2],dt,y[0],tbegin,tend,
 332                        oenv,TRUE,eFitFn,fitparms,fix);
 333   for(i=0; i<nx; i++)
 334     y[3][i] = fit_function(eFitFn,fitparms,y[0][i]);
 335
 336   if (epsRF == 0) {
 337     /* This means infinity! */
 338     lambda = 0;
 339     rffac  = 1;
 340   }
 341   else {
 342     lambda = (eps0 - 1.0)/(2*epsRF - 1.0);
 343     rffac  = (2*epsRF+eps0)/(2*epsRF+1);
 344   }
 345   printf("DATA INTEGRAL: %5.1f, tauD(old) = %5.1f ps, "
 346           "tau_slope = %5.1f, tau_slope,D = %5.1f ps\n",
 347           integral,integral*rffac,fitparms[0],fitparms[0]*rffac);
 348
 349   printf("tau_D from tau1 = %8.3g , eps(Infty) = %8.3f\n",
 350           fitparms[0]*(1 + fitparms[1]*lambda),
 351           1 + ((1 - fitparms[1])*(eps0 - 1))/(1 + fitparms[1]*lambda));
 352
 353   fitintegral=numerical_deriv(nx,y[0],y[1],y[3],y[4],y[5],tendInt,nsmooth);
 354   printf("FIT INTEGRAL (tau_M): %5.1f, tau_D = %5.1f\n",
 355           fitintegral,fitintegral*rffac);
 356
 357   /* Now we have the negative gradient of <Phi(0) Phi(t)> */
 358   write_xvg(opt2fn("-d",NFILE,fnm),"Data",nx-1,6,y,legend,oenv);
 359
 360   /* Do FFT and analysis */
 361   do_four(opt2fn("-o",NFILE,fnm),opt2fn("-c",NFILE,fnm),
 362           nx-1,y[0],y[5],eps0,epsRF,oenv);
 363
 364   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 365   do_view(oenv,opt2fn("-c",NFILE,fnm),NULL);
 366   do_view(oenv,opt2fn("-d",NFILE,fnm),"-nxy");
 367
 368   thanx(stderr);
 369
 370   return 0;
 371 }