src/gromacs/gmxana/gmx_dielectric.cpp

   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  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source 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 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdio>
  42 #include <cstdlib>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/commandline/viewit.h"
  48 #include "gromacs/correlationfunctions/expfit.h"
  49 #include "gromacs/correlationfunctions/integrate.h"
  50 #include "gromacs/fft/fft.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gromacs/gmxana/gmx_ana.h"
  53 #include "gromacs/gmxana/gstat.h"
  54 #include "gromacs/math/gmxcomplex.h"
  55 #include "gromacs/math/utilities.h"
  56 #include "gromacs/utility/arraysize.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/futil.h"
  59 #include "gromacs/utility/pleasecite.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 /* Determines at which point in the array the fit should start */
  63 static int calc_nbegin(int nx, real x[], real tbegin)
  64 {
  65     int nbegin;
  66
  67     /* Assume input x is sorted */
  68     for (nbegin = 0; (nbegin < nx) && (x[nbegin] <= tbegin); nbegin++)
  69     {
  70         // Deliberately left empty
  71     }
  72     if ((nbegin == nx) || (nbegin == 0))
  73     {
  74         gmx_fatal(FARGS, "Begin time %f not in x-domain [%f through %f]\n", tbegin, x[0], x[nx - 1]);
  75     }
  76
  77     /* Take the one closest to tbegin */
  78     if (std::abs(x[nbegin] - tbegin) > std::abs(x[nbegin - 1] - tbegin))
  79     {
  80         nbegin--;
  81     }
  82
  83     printf("nbegin = %d, x[nbegin] = %g, tbegin = %g\n", nbegin, x[nbegin], tbegin);
  84
  85     return nbegin;
  86 }
  87
  88 static real
  89 numerical_deriv(int nx, real x[], const real y[], const real fity[], real combined[], real dy[], real tendInt, int nsmooth)
  90 {
  91     FILE* tmpfp;
  92     int   i, nbegin, i0, i1;
  93     real  fac, fx, fy, integralSmth;
  94
  95     nbegin = calc_nbegin(nx, x, tendInt);
  96     if (nsmooth == 0)
  97     {
  98         for (i = 0; (i < nbegin); i++)
  99         {
 100             combined[i] = y[i];
 101         }
 102         fac = y[nbegin] / fity[nbegin];
 103         printf("scaling fitted curve by %g\n", fac);
 104         for (i = nbegin; (i < nx); i++)
 105         {
 106             combined[i] = fity[i] * fac;
 107         }
 108     }
 109     else
 110     {
 111         i0 = std::max(0, nbegin);
 112         i1 = std::min(nx - 1, nbegin + nsmooth);
 113         printf("Making smooth transition from %d through %d\n", i0, i1);
 114         for (i = 0; (i < i0); i++)
 115         {
 116             combined[i] = y[i];
 117         }
 118         for (i = i0; (i <= i1); i++)
 119         {
 120             fx = static_cast<real>(i1 - i) / (i1 - i0);
 121             fy = static_cast<real>(i - i0) / (i1 - i0);
 122             if (debug)
 123             {
 124                 fprintf(debug, "x: %g factors for smoothing: %10g %10g\n", x[i], fx, fy);
 125             }
 126             combined[i] = fx * y[i] + fy * fity[i];
 127         }
 128         for (i = i1 + 1; (i < nx); i++)
 129         {
 130             combined[i] = fity[i];
 131         }
 132     }
 133
 134     tmpfp        = gmx_ffopen("integral_smth.xvg", "w");
 135     integralSmth = print_and_integrate(tmpfp, nx, x[1] - x[0], combined, nullptr, 1);
 136     printf("SMOOTH integral = %10.5e\n", integralSmth);
 137
 138     dy[0] = (combined[1] - combined[0]) / (x[1] - x[0]);
 139     for (i = 1; (i < nx - 1); i++)
 140     {
 141         dy[i] = (combined[i + 1] - combined[i - 1]) / (x[i + 1] - x[i - 1]);
 142     }
 143     dy[nx - 1] = (combined[nx - 1] - combined[nx - 2]) / (x[nx - 1] - x[nx - 2]);
 144
 145     for (i = 0; (i < nx); i++)
 146     {
 147         dy[i] *= -1;
 148     }
 149
 150     return integralSmth;
 151 }
 152
 153 static void do_four(const char*             fn,
 154                     const char*             cn,
 155                     int                     nx,
 156                     const real              x[],
 157                     const real              dy[],
 158                     real                    eps0,
 159                     real                    epsRF,
 160                     const gmx_output_env_t* oenv)
 161 {
 162     FILE *     fp, *cp;
 163     t_complex *tmp, gw, hw, kw;
 164     int        i, nnx, nxsav;
 165     real       fac, nu, dt, maxeps, numax;
 166     gmx_fft_t  fft;
 167     int        fftcode;
 168
 169     nxsav = nx;
 170     /*while ((dy[nx-1] == 0.0) && (nx > 0))
 171        nx--;*/
 172     if (nx == 0)
 173     {
 174         gmx_fatal(FARGS, "Empty dataset in %s, line %d!", __FILE__, __LINE__);
 175     }
 176
 177     nnx = 1;
 178     while (nnx < 2 * nx)
 179     {
 180         nnx *= 2;
 181     }
 182
 183     snew(tmp, 2 * nnx);
 184     printf("Doing FFT of %d points\n", nnx);
 185     for (i = 0; (i < nx); i++)
 186     {
 187         tmp[i].re = dy[i];
 188     }
 189     if ((fftcode = gmx_fft_init_1d_real(&fft, nnx, GMX_FFT_FLAG_NONE)) != 0)
 190     {
 191         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 192     }
 193     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_COMPLEX_TO_REAL, tmp, tmp)) != 0)
 194     {
 195         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 196     }
 197     gmx_fft_destroy(fft);
 198     dt = x[1] - x[0];
 199     if (epsRF == 0)
 200     {
 201         fac = (eps0 - 1) / tmp[0].re;
 202     }
 203     else
 204     {
 205         fac = ((eps0 - 1) / (2 * epsRF + eps0)) / tmp[0].re;
 206     }
 207     fp     = xvgropen(fn, "Epsilon(\\8w\\4)", "Freq. (GHz)", "eps", oenv);
 208     cp     = xvgropen(cn, "Cole-Cole plot", "Eps'", "Eps''", oenv);
 209     maxeps = 0;
 210     numax  = 0;
 211     for (i = 0; (i < nxsav); i++)
 212     {
 213         if (epsRF == 0)
 214         {
 215             kw.re = 1 + fac * tmp[i].re;
 216             kw.im = 1 + fac * tmp[i].im;
 217         }
 218         else
 219         {
 220             gw = rcmul(fac, tmp[i]);
 221             hw = rcmul(2 * epsRF, gw);
 222             hw.re += 1.0;
 223             gw.re = 1.0 - gw.re;
 224             gw.im = -gw.im;
 225             kw    = cdiv(hw, gw);
 226         }
 227         kw.im *= -1;
 228
 229         nu = (i + 1) * 1000.0 / (nnx * dt);
 230         if (kw.im > maxeps)
 231         {
 232             maxeps = kw.im;
 233             numax  = nu;
 234         }
 235
 236         fprintf(fp, "%10.5e  %10.5e  %10.5e\n", nu, kw.re, kw.im);
 237         fprintf(cp, "%10.5e  %10.5e\n", kw.re, kw.im);
 238     }
 239     printf("MAXEPS = %10.5e at frequency %10.5e GHz (tauD = %8.1f ps)\n", maxeps, numax,
 240            1000 / (2 * M_PI * numax));
 241     xvgrclose(fp);
 242     xvgrclose(cp);
 243     sfree(tmp);
 244 }
 245
 246 int gmx_dielectric(int argc, char* argv[])
 247 {
 248     const char* desc[] = {
 249         "[THISMODULE] calculates frequency dependent dielectric constants",
 250         "from the autocorrelation function of the total dipole moment in",
 251         "your simulation. This ACF can be generated by [gmx-dipoles].",
 252         "The functional forms of the available functions are:",
 253         "",
 254         " * One parameter:    y = [EXP]-a[SUB]1[sub] x[exp],",
 255         " * Two parameters:   y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp],",
 256         " * Three parameters: y = a[SUB]2[sub] [EXP]-a[SUB]1[sub] x[exp] + (1 - ",
 257         "   a[SUB]2[sub]) [EXP]-a[SUB]3[sub] x[exp].",
 258         "",
 259         "Start values for the fit procedure can be given on the command line.",
 260         "It is also possible to fix parameters at their start value, use [TT]-fix[tt]",
 261         "with the number of the parameter you want to fix.",
 262         "[PAR]",
 263         "Three output files are generated, the first contains the ACF,",
 264         "an exponential fit to it with 1, 2 or 3 parameters, and the",
 265         "numerical derivative of the combination data/fit.",
 266         "The second file contains the real and imaginary parts of the",
 267         "frequency-dependent dielectric constant, the last gives a plot",
 268         "known as the Cole-Cole plot, in which the imaginary",
 269         "component is plotted as a function of the real component.",
 270         "For a pure exponential relaxation (Debye relaxation) the latter",
 271         "plot should be one half of a circle."
 272     };
 273     t_filenm fnm[] = { { efXVG, "-f", "dipcorr", ffREAD },
 274                        { efXVG, "-d", "deriv", ffWRITE },
 275                        { efXVG, "-o", "epsw", ffWRITE },
 276                        { efXVG, "-c", "cole", ffWRITE } };
 277 #define NFILE asize(fnm)
 278     gmx_output_env_t* oenv;
 279     int               i, j, nx, ny, nxtail, eFitFn, nfitparm;
 280     real              dt, integral, fitintegral, fac, rffac;
 281     double*           fitparms;
 282     double**          yd;
 283     real**            y;
 284     const char*       legend[] = { "Correlation", "Std. Dev.", "Fit", "Combined", "Derivative" };
 285     static int        fix = 0, bX = 1, nsmooth = 3;
 286     static real       tendInt = 5.0, tbegin = 5.0, tend = 500.0;
 287     static real       A = 0.5, tau1 = 10.0, tau2 = 1.0, eps0 = 80, epsRF = 78.5, tail = 500.0;
 288     real              lambda;
 289     t_pargs           pa[] = {
 290         { "-x1",
 291           FALSE,
 292           etBOOL,
 293           { &bX },
 294           "use first column as [IT]x[it]-axis rather than first data set" },
 295         { "-eint",
 296           FALSE,
 297           etREAL,
 298           { &tendInt },
 299           "Time to end the integration of the data and start to use the fit" },
 300         { "-bfit", FALSE, etREAL, { &tbegin }, "Begin time of fit" },
 301         { "-efit", FALSE, etREAL, { &tend }, "End time of fit" },
 302         { "-tail", FALSE, etREAL, { &tail }, "Length of function including data and tail from fit" },
 303         { "-A", FALSE, etREAL, { &A }, "Start value for fit parameter A" },
 304         { "-tau1", FALSE, etREAL, { &tau1 }, "Start value for fit parameter [GRK]tau[grk]1" },
 305         { "-tau2", FALSE, etREAL, { &tau2 }, "Start value for fit parameter [GRK]tau[grk]2" },
 306         { "-eps0", FALSE, etREAL, { &eps0 }, "[GRK]epsilon[grk]0 of your liquid" },
 307         { "-epsRF",
 308           FALSE,
 309           etREAL,
 310           { &epsRF },
 311           "[GRK]epsilon[grk] of the reaction field used in your simulation. A value of 0 means "
 312           "infinity." },
 313         { "-fix",
 314           FALSE,
 315           etINT,
 316           { &fix },
 317           "Fix parameters at their start values, A (2), tau1 (1), or tau2 (4)" },
 318         { "-ffn", FALSE, etENUM, { s_ffn }, "Fit function" },
 319         { "-nsmooth", FALSE, etINT, { &nsmooth }, "Number of points for smoothing" }
 320     };
 321
 322     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW, NFILE, fnm, asize(pa), pa,
 323                            asize(desc), desc, 0, nullptr, &oenv))
 324     {
 325         return 0;
 326     }
 327     please_cite(stdout, "Spoel98a");
 328     printf("WARNING: non-polarizable models can never yield an infinite\n"
 329            "dielectric constant that is different from 1. This is incorrect\n"
 330            "in the reference given above (Spoel98a).\n\n");
 331
 332
 333     nx     = read_xvg(opt2fn("-f", NFILE, fnm), &yd, &ny);
 334     dt     = yd[0][1] - yd[0][0];
 335     nxtail = std::min(static_cast<int>(tail / dt), nx);
 336
 337     printf("Read data set containing %d colums and %d rows\n", ny, nx);
 338     printf("Assuming (from data) that timestep is %g, nxtail = %d\n", dt, nxtail);
 339     snew(y, 6);
 340     for (i = 0; (i < ny); i++)
 341     {
 342         snew(y[i], std::max(nx, nxtail));
 343     }
 344     for (i = 0; (i < nx); i++)
 345     {
 346         y[0][i] = yd[0][i];
 347         for (j = 1; (j < ny); j++)
 348         {
 349             y[j][i] = yd[j][i];
 350         }
 351     }
 352     if (nxtail > nx)
 353     {
 354         for (i = nx; (i < nxtail); i++)
 355         {
 356             y[0][i] = dt * i + y[0][0];
 357             for (j = 1; (j < ny); j++)
 358             {
 359                 y[j][i] = 0.0;
 360             }
 361         }
 362         nx = nxtail;
 363     }
 364
 365
 366     /* We have read a file WITHOUT standard deviations, so we make our own... */
 367     if (ny == 2)
 368     {
 369         printf("Creating standard deviation numbers ...\n");
 370         snew(y[2], nx);
 371
 372         fac = 1.0 / nx;
 373         for (i = 0; (i < nx); i++)
 374         {
 375             y[2][i] = fac;
 376         }
 377     }
 378
 379     eFitFn   = sffn2effn(s_ffn);
 380     nfitparm = effnNparams(eFitFn);
 381     snew(fitparms, 4);
 382     fitparms[0] = tau1;
 383     if (nfitparm > 1)
 384     {
 385         fitparms[1] = A;
 386     }
 387     if (nfitparm > 2)
 388     {
 389         fitparms[2] = tau2;
 390     }
 391
 392
 393     snew(y[3], nx);
 394     snew(y[4], nx);
 395     snew(y[5], nx);
 396
 397     integral = print_and_integrate(nullptr, calc_nbegin(nx, y[0], tbegin), dt, y[1], nullptr, 1);
 398     integral += do_lmfit(nx, y[1], y[2], dt, y[0], tbegin, tend, oenv, TRUE, eFitFn, fitparms, fix, nullptr);
 399     for (i = 0; i < nx; i++)
 400     {
 401         y[3][i] = fit_function(eFitFn, fitparms, y[0][i]);
 402     }
 403
 404     if (epsRF == 0)
 405     {
 406         /* This means infinity! */
 407         lambda = 0;
 408         rffac  = 1;
 409     }
 410     else
 411     {
 412         lambda = (eps0 - 1.0) / (2 * epsRF - 1.0);
 413         rffac  = (2 * epsRF + eps0) / (2 * epsRF + 1);
 414     }
 415     printf("DATA INTEGRAL: %5.1f, tauD(old) = %5.1f ps, "
 416            "tau_slope = %5.1f, tau_slope,D = %5.1f ps\n",
 417            integral, integral * rffac, fitparms[0], fitparms[0] * rffac);
 418
 419     printf("tau_D from tau1 = %8.3g , eps(Infty) = %8.3f\n", fitparms[0] * (1 + fitparms[1] * lambda),
 420            1 + ((1 - fitparms[1]) * (eps0 - 1)) / (1 + fitparms[1] * lambda));
 421
 422     fitintegral = numerical_deriv(nx, y[0], y[1], y[3], y[4], y[5], tendInt, nsmooth);
 423     printf("FIT INTEGRAL (tau_M): %5.1f, tau_D = %5.1f\n", fitintegral, fitintegral * rffac);
 424
 425     /* Now we have the negative gradient of <Phi(0) Phi(t)> */
 426     write_xvg(opt2fn("-d", NFILE, fnm), "Data", nx - 1, 6, y, legend, oenv);
 427
 428     /* Do FFT and analysis */
 429     do_four(opt2fn("-o", NFILE, fnm), opt2fn("-c", NFILE, fnm), nx - 1, y[0], y[5], eps0, epsRF, oenv);
 430
 431     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 432     do_view(oenv, opt2fn("-c", NFILE, fnm), nullptr);
 433     do_view(oenv, opt2fn("-d", NFILE, fnm), "-nxy");
 434
 435     return 0;
 436 }