src/gromacs/correlationfunctions/expfit.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-2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal
  38  * \file
  39  * \brief
  40  * Implements routine for fitting a data set to a curve
  41  *
  42  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  43  * \ingroup module_correlationfunctions
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "expfit.h"
  48
  49 #include <cstring>
  50
  51 #include <algorithm>
  52
  53 #include "gromacs/correlationfunctions/integrate.h"
  54 #include "gromacs/fileio/xvgr.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/futil.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/real.h"
  60 #include "gromacs/utility/smalloc.h"
  61
  62 #include "gmx_lmcurve.h"
  63
  64 /*! \brief Number of parameters for each fitting function */
  65 static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 5, 7, 9, 2, 4, 3, 6 };
  66
  67 /* +2 becuase parse_common_args wants leading and concluding NULL.
  68  * We only allow exponential functions as choices on the command line,
  69  * hence there are many more NULL field (which have to be at the end of
  70  * the array).
  71  */
  72 const char* s_ffn[effnNR + 2] = { nullptr, "none",  "exp",   "aexp",  "exp_exp", "exp5", "exp7",
  73                                   "exp9",  nullptr, nullptr, nullptr, nullptr,   nullptr };
  74
  75 // clang-format off
  76 // needed because clang-format wants to break the lines below
  77 /*! \brief Long description for each fitting function type */
  78 static const char* longs_ffn[effnNR]
  79         = { "no fit",
  80             "y = exp(-x/|a0|)",
  81             "y = a1 exp(-x/|a0|)",
  82             "y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0",
  83             "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4, a3 >= a1",
  84             "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1",
  85             "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1",
  86             "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)",
  87             "y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)",
  88             "y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)",
  89             "y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)" };
  90 // clang-format on
  91
  92 int effnNparams(int effn)
  93 {
  94     if ((0 <= effn) && (effn < effnNR))
  95     {
  96         return nfp_ffn[effn];
  97     }
  98     else
  99     {
 100         return -1;
 101     }
 102 }
 103
 104 const char* effnDescription(int effn)
 105 {
 106     if ((0 <= effn) && (effn < effnNR))
 107     {
 108         return longs_ffn[effn];
 109     }
 110     else
 111     {
 112         return nullptr;
 113     }
 114 }
 115
 116 int sffn2effn(const char** sffn)
 117 {
 118     int eFitFn, i;
 119
 120     eFitFn = 0;
 121     for (i = 0; i < effnNR; i++)
 122     {
 123         if (sffn[i + 1] && strcmp(sffn[0], sffn[i + 1]) == 0)
 124         {
 125             eFitFn = i;
 126         }
 127     }
 128
 129     return eFitFn;
 130 }
 131
 132 /*! \brief Compute exponential function A exp(-x/tau) */
 133 static double myexp(double x, double A, double tau)
 134 {
 135     if ((A == 0) || (tau == 0))
 136     {
 137         return 0;
 138     }
 139     return A * exp(-x / tau);
 140 }
 141
 142 /*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
 143 static double lmc_erffit(double x, const double* a)
 144 {
 145     double erfarg;
 146     double myerf;
 147
 148     if (a[3] != 0)
 149     {
 150         erfarg = (x - a[2]) / (a[3] * a[3]);
 151         myerf  = std::erf(erfarg);
 152     }
 153     else
 154     {
 155         /* If a[3] == 0, a[3]^2 = 0 and the erfarg becomes +/- infinity */
 156         if (x < a[2])
 157         {
 158             myerf = -1;
 159         }
 160         else
 161         {
 162             myerf = 1;
 163         }
 164     }
 165     return 0.5 * ((a[0] + a[1]) - (a[0] - a[1]) * myerf);
 166 }
 167
 168 /*! \brief Exponent function that prevents overflow */
 169 static double safe_exp(double x)
 170 {
 171     double exp_max = 200;
 172     double exp_min = -exp_max;
 173     if (x <= exp_min)
 174     {
 175         return exp(exp_min);
 176     }
 177     else if (x >= exp_max)
 178     {
 179         return exp(exp_max);
 180     }
 181     else
 182     {
 183         return exp(x);
 184     }
 185 }
 186
 187 /*! \brief Exponent minus 1 function that prevents overflow */
 188 static double safe_expm1(double x)
 189 {
 190     double exp_max = 200;
 191     double exp_min = -exp_max;
 192     if (x <= exp_min)
 193     {
 194         return -1;
 195     }
 196     else if (x >= exp_max)
 197     {
 198         return exp(exp_max);
 199     }
 200     else
 201     {
 202         return std::expm1(x);
 203     }
 204 }
 205
 206 /*! \brief Compute y = exp(-x/|a0|) */
 207 static double lmc_exp_one_parm(double x, const double* a)
 208 {
 209     return safe_exp(-x / fabs(a[0]));
 210 }
 211
 212 /*! \brief Compute y = a1 exp(-x/|a0|) */
 213 static double lmc_exp_two_parm(double x, const double* a)
 214 {
 215     return a[1] * safe_exp(-x / fabs(a[0]));
 216 }
 217
 218 /*! \brief Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|) */
 219 static double lmc_exp_exp(double x, const double* a)
 220 {
 221     double e1, e2;
 222
 223     e1 = safe_exp(-x / fabs(a[0]));
 224     e2 = safe_exp(-x / (fabs(a[0]) + fabs(a[2])));
 225     return a[1] * e1 + (1 - a[1]) * e2;
 226 }
 227
 228 /*! \brief Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 */
 229 static double lmc_exp_5_parm(double x, const double* a)
 230 {
 231     double e1, e2;
 232
 233     e1 = safe_exp(-x / fabs(a[1]));
 234     e2 = safe_exp(-x / (fabs(a[1]) + fabs(a[3])));
 235     return a[0] * e1 + a[2] * e2 + a[4];
 236 }
 237
 238 /*! \brief Compute 7 parameter exponential function value.
 239  *
 240  * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
 241  * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6
 242  */
 243 static double lmc_exp_7_parm(double x, const double* a)
 244 {
 245     double e1, e2, e3;
 246     double fa1, fa3, fa5;
 247
 248     fa1 = fabs(a[1]);
 249     fa3 = fa1 + fabs(a[3]);
 250     fa5 = fa3 + fabs(a[5]);
 251     e1  = safe_exp(-x / fa1);
 252     e2  = safe_exp(-x / fa3);
 253     e3  = safe_exp(-x / fa5);
 254     return a[0] * e1 + a[2] * e2 + a[4] * e3 + a[6];
 255 }
 256
 257 /*! \brief Compute 9 parameter exponential function value.
 258  *
 259  * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
 260  * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8
 261  */
 262 static double lmc_exp_9_parm(double x, const double* a)
 263 {
 264     double e1, e2, e3, e4;
 265     double fa1, fa3, fa5, fa7;
 266
 267     fa1 = fabs(a[1]);
 268     fa3 = fa1 + fabs(a[3]);
 269     fa5 = fa3 + fabs(a[5]);
 270     fa7 = fa5 + fabs(a[7]);
 271
 272     e1 = safe_exp(-x / fa1);
 273     e2 = safe_exp(-x / fa3);
 274     e3 = safe_exp(-x / fa5);
 275     e4 = safe_exp(-x / fa7);
 276     return a[0] * e1 + a[2] * e2 + a[4] * e3 + a[6] * e4 + a[8];
 277 }
 278
 279 /*! \brief Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|) */
 280 static double lmc_pres_6_parm(double x, const double* a)
 281 {
 282     double term1, term2, term3;
 283     double pow_max = 10;
 284
 285     term3 = 0;
 286     if ((a[4] != 0) && (a[0] != 0))
 287     {
 288         double power = std::min(fabs(a[5]), pow_max);
 289         term3        = a[0] * safe_exp(-pow((x / fabs(a[4])), power));
 290     }
 291
 292     term1 = 1 - a[0];
 293     term2 = 0;
 294     if ((term1 != 0) && (a[2] != 0))
 295     {
 296         double power = std::min(fabs(a[3]), pow_max);
 297         term2        = safe_exp(-pow((x / fabs(a[2])), power)) * cos(x * fabs(a[1]));
 298     }
 299
 300     return term1 * term2 + term3;
 301 }
 302
 303 /*! \brief Compute vac function */
 304 static double lmc_vac_2_parm(double x, const double* a)
 305 {
 306     /* Fit to function
 307      *
 308      * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 309      *
 310      *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 311      *
 312      *    v = x/(2 a0)
 313      *    w = sqrt(1 - a1)
 314      *
 315      *    For tranverse current autocorrelation functions:
 316      *       a0 = tau
 317      *       a1 = 4 tau (eta/rho) k^2
 318      *
 319      */
 320
 321     double y, v, det, omega, wv, em, ec, es;
 322     double wv_max = 100;
 323
 324     v   = x / (2 * fabs(a[0]));
 325     det = 1 - a[1];
 326     em  = safe_exp(-v);
 327     if (det != 0)
 328     {
 329         omega = sqrt(fabs(det));
 330         wv    = std::min(omega * v, wv_max);
 331
 332         if (det > 0)
 333         {
 334             ec = em * 0.5 * (safe_exp(wv) + safe_exp(-wv));
 335             es = em * 0.5 * (safe_exp(wv) - safe_exp(-wv)) / omega;
 336         }
 337         else
 338         {
 339             ec = em * cos(wv);
 340             es = em * sin(wv) / omega;
 341         }
 342         y = ec + es;
 343     }
 344     else
 345     {
 346         y = (1 + v) * em;
 347     }
 348     return y;
 349 }
 350
 351 /*! \brief Compute error estimate */
 352 static double lmc_errest_3_parm(double x, const double* a)
 353 {
 354     double e1, e2, v1, v2;
 355     double fa0 = fabs(a[0]);
 356     double fa1;
 357     double fa2 = fa0 + fabs(a[2]);
 358
 359     if (a[0] != 0)
 360     {
 361         e1 = safe_expm1(-x / fa0);
 362     }
 363     else
 364     {
 365         e1 = 0;
 366     }
 367     if (a[2] != 0)
 368     {
 369         e2 = safe_expm1(-x / fa2);
 370     }
 371     else
 372     {
 373         e2 = 0;
 374     }
 375
 376     if (x > 0)
 377     {
 378         v1 = 2 * fa0 * (e1 * fa0 / x + 1);
 379         v2 = 2 * fa2 * (e2 * fa2 / x + 1);
 380         /* We need 0 <= a1 <= 1 */
 381         fa1 = std::min(1.0, std::max(0.0, a[1]));
 382
 383         return fa1 * v1 + (1 - fa1) * v2;
 384     }
 385     else
 386     {
 387         return 0;
 388     }
 389 }
 390
 391 /*! \brief array of fitting functions corresponding to the pre-defined types */
 392 t_lmcurve lmcurves[effnNR + 1] = { lmc_exp_one_parm,  lmc_exp_one_parm, lmc_exp_two_parm,
 393                                    lmc_exp_exp,       lmc_exp_5_parm,   lmc_exp_7_parm,
 394                                    lmc_exp_9_parm,    lmc_vac_2_parm,   lmc_erffit,
 395                                    lmc_errest_3_parm, lmc_pres_6_parm };
 396
 397 double fit_function(const int eFitFn, const double parm[], const double x)
 398 {
 399     if ((eFitFn < 0) || (eFitFn >= effnNR))
 400     {
 401         fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", eFitFn, effnNR - 1);
 402         return 0.0;
 403     }
 404     return lmcurves[eFitFn](x, parm);
 405 }
 406
 407 /*! \brief Ensure the fitting parameters are well-behaved.
 408  *
 409  * In order to make sure that fitting behaves according to the
 410  * constraint that time constants are positive and increasing
 411  * we enforce this here by setting all time constants to their
 412  * absolute value and by adding e.g. |a_0| to |a_2|. This is
 413  * done in conjunction with the fitting functions themselves.
 414  * When there are multiple time constants we make sure that
 415  * the successive time constants are at least double the
 416  * previous ones and during fitting we enforce the they remain larger.
 417  * It may very well help the convergence of the fitting routine.
 418  */
 419 static void initiate_fit_params(int eFitFn, double params[])
 420 {
 421     int i, nparm;
 422
 423     nparm = effnNparams(eFitFn);
 424
 425     switch (eFitFn)
 426     {
 427         case effnVAC: GMX_ASSERT(params[0] >= 0, "parameters should be >= 0"); break;
 428         case effnEXP1:
 429         case effnEXP2:
 430         case effnEXPEXP:
 431             GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
 432             if (nparm > 2)
 433             {
 434                 GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
 435                 /* In order to maintain params[2] >= params[0] in the final
 436                  * result, we fit the difference between the two, that is
 437                  * params[2]-params[0] and in the add add in params[0]
 438                  * again.
 439                  */
 440                 params[2] = std::max(fabs(params[2]) - params[0], params[0]);
 441             }
 442             break;
 443         case effnEXP5:
 444         case effnEXP7:
 445         case effnEXP9:
 446             GMX_ASSERT(params[1] >= 0, "parameters should be >= 0");
 447             params[1] = fabs(params[1]);
 448             if (nparm > 3)
 449             {
 450                 GMX_ASSERT(params[3] >= 0, "parameters should be >= 0");
 451                 /* See comment under effnEXPEXP */
 452                 params[3] = std::max(fabs(params[3]) - params[1], params[1]);
 453                 if (nparm > 5)
 454                 {
 455                     GMX_ASSERT(params[5] >= 0, "parameters should be >= 0");
 456                     /* See comment under effnEXPEXP */
 457                     params[5] = std::max(fabs(params[5]) - params[3], params[3]);
 458                     if (nparm > 7)
 459                     {
 460                         GMX_ASSERT(params[7] >= 0, "parameters should be >= 0");
 461                         /* See comment under effnEXPEXP */
 462                         params[7] = std::max(fabs(params[7]) - params[5], params[5]);
 463                     }
 464                 }
 465             }
 466             break;
 467         case effnERREST:
 468             GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
 469             GMX_ASSERT(params[1] >= 0 && params[1] <= 1, "parameter 1 should in 0 .. 1");
 470             GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
 471             /* See comment under effnEXPEXP */
 472             params[2] = fabs(params[2]) - params[0];
 473             break;
 474         case effnPRES:
 475             for (i = 1; (i < nparm); i++)
 476             {
 477                 GMX_ASSERT(params[i] >= 0, "parameters should be >= 0");
 478             }
 479             break;
 480         default: break;
 481     }
 482 }
 483
 484 /*! \brief Process the fitting parameters to get output parameters.
 485  *
 486  * See comment at the previous function.
 487  */
 488 static void extract_fit_params(int eFitFn, double params[])
 489 {
 490     int i, nparm;
 491
 492     nparm = effnNparams(eFitFn);
 493
 494     switch (eFitFn)
 495     {
 496         case effnVAC: params[0] = fabs(params[0]); break;
 497         case effnEXP1:
 498         case effnEXP2:
 499         case effnEXPEXP:
 500             params[0] = fabs(params[0]);
 501             if (nparm > 2)
 502             {
 503                 /* Back conversion of parameters from the fitted difference
 504                  * to the absolute value.
 505                  */
 506                 params[2] = fabs(params[2]) + params[0];
 507             }
 508             break;
 509         case effnEXP5:
 510         case effnEXP7:
 511         case effnEXP9:
 512             params[1] = fabs(params[1]);
 513             if (nparm > 3)
 514             {
 515                 /* See comment under effnEXPEXP */
 516                 params[3] = fabs(params[3]) + params[1];
 517                 if (nparm > 5)
 518                 {
 519                     /* See comment under effnEXPEXP */
 520                     params[5] = fabs(params[5]) + params[3];
 521                     if (nparm > 7)
 522                     {
 523                         /* See comment under effnEXPEXP */
 524                         params[7] = fabs(params[7]) + params[5];
 525                     }
 526                 }
 527             }
 528             break;
 529         case effnERREST:
 530             params[0] = fabs(params[0]);
 531             if (params[1] < 0)
 532             {
 533                 params[1] = 0;
 534             }
 535             else if (params[1] > 1)
 536             {
 537                 params[1] = 1;
 538             }
 539             /* See comment under effnEXPEXP */
 540             params[2] = params[0] + fabs(params[2]);
 541             break;
 542         case effnPRES:
 543             for (i = 1; (i < nparm); i++)
 544             {
 545                 params[i] = fabs(params[i]);
 546             }
 547             break;
 548         default: break;
 549     }
 550 }
 551
 552 /*! \brief Print chi-squared value and the parameters */
 553 static void print_chi2_params(FILE*        fp,
 554                               const int    eFitFn,
 555                               const double fitparms[],
 556                               const char*  label,
 557                               const int    nfitpnts,
 558                               const double x[],
 559                               const double y[])
 560 {
 561     int    i;
 562     double chi2 = 0;
 563
 564     for (i = 0; (i < nfitpnts); i++)
 565     {
 566         double yfit = lmcurves[eFitFn](x[i], fitparms);
 567         chi2 += gmx::square(y[i] - yfit);
 568     }
 569     fprintf(fp,
 570             "There are %d data points, %d parameters, %s chi2 = %g\nparams:",
 571             nfitpnts,
 572             effnNparams(eFitFn),
 573             label,
 574             chi2);
 575     for (i = 0; (i < effnNparams(eFitFn)); i++)
 576     {
 577         fprintf(fp, "  %10g", fitparms[i]);
 578     }
 579     fprintf(fp, "\n");
 580 }
 581
 582 real do_lmfit(int                     ndata,
 583               const real              c1[],
 584               real                    sig[],
 585               real                    dt,
 586               const real*             x0,
 587               real                    begintimefit,
 588               real                    endtimefit,
 589               const gmx_output_env_t* oenv,
 590               gmx_bool                bVerbose,
 591               int                     eFitFn,
 592               double                  fitparms[],
 593               int                     fix,
 594               const char*             fn_fitted)
 595 {
 596     FILE*   fp;
 597     int     i, j, nfitpnts;
 598     double  integral, ttt;
 599     double *x, *y, *dy;
 600
 601     if (0 != fix)
 602     {
 603         fprintf(stderr, "Using fixed parameters in curve fitting is temporarily not working.\n");
 604     }
 605     if (debug)
 606     {
 607         fprintf(debug, "There are %d points to fit %d vars!\n", ndata, effnNparams(eFitFn));
 608         fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n", eFitFn, begintimefit, endtimefit, dt);
 609     }
 610
 611     snew(x, ndata);
 612     snew(y, ndata);
 613     snew(dy, ndata);
 614     j = 0;
 615     for (i = 0; i < ndata; i++)
 616     {
 617         ttt = x0 ? x0[i] : dt * i;
 618         if (ttt >= begintimefit && ttt <= endtimefit)
 619         {
 620             x[j] = ttt;
 621             y[j] = c1[i];
 622             if (nullptr == sig)
 623             {
 624                 // No weighting if all values are divided by 1.
 625                 dy[j] = 1;
 626             }
 627             else
 628             {
 629                 dy[j] = std::max(1.0e-7, static_cast<double>(sig[i]));
 630             }
 631             if (debug)
 632             {
 633                 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy=%g, ttt=%g\n", j, i, x[j], y[j], dy[j], ttt);
 634             }
 635             j++;
 636         }
 637     }
 638     nfitpnts = j;
 639     integral = 0;
 640     if (nfitpnts < effnNparams(eFitFn))
 641     {
 642         fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!\n", nfitpnts, dt);
 643     }
 644     else
 645     {
 646         gmx_bool bSuccess;
 647
 648         if (bVerbose)
 649         {
 650             print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y);
 651         }
 652         initiate_fit_params(eFitFn, fitparms);
 653
 654         bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix);
 655         extract_fit_params(eFitFn, fitparms);
 656
 657         if (!bSuccess)
 658         {
 659             fprintf(stderr, "Fit failed!\n");
 660         }
 661         else
 662         {
 663             if (bVerbose)
 664             {
 665                 print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y);
 666             }
 667             switch (eFitFn)
 668             {
 669                 case effnEXP1: integral = fitparms[0] * myexp(begintimefit, 1, fitparms[0]); break;
 670                 case effnEXP2:
 671                     integral = fitparms[0] * myexp(begintimefit, fitparms[1], fitparms[0]);
 672                     break;
 673                 case effnEXPEXP:
 674                     integral = (fitparms[0] * myexp(begintimefit, fitparms[1], fitparms[0])
 675                                 + fitparms[2] * myexp(begintimefit, 1 - fitparms[1], fitparms[2]));
 676                     break;
 677                 case effnEXP5:
 678                 case effnEXP7:
 679                 case effnEXP9:
 680                     integral = 0;
 681                     for (i = 0; (i < (effnNparams(eFitFn) - 1) / 2); i++)
 682                     {
 683                         integral += fitparms[2 * i]
 684                                     * myexp(begintimefit, fitparms[2 * i + 1], fitparms[2 * i]);
 685                     }
 686                     break;
 687                 default:
 688                     /* Do numerical integration */
 689                     integral = 0;
 690                     for (i = 0; (i < nfitpnts - 1); i++)
 691                     {
 692                         double y0 = lmcurves[eFitFn](x[i], fitparms);
 693                         double y1 = lmcurves[eFitFn](x[i + 1], fitparms);
 694                         integral += (x[i + 1] - x[i]) * (y1 + y0) * 0.5;
 695                     }
 696                     break;
 697             }
 698
 699             if (bVerbose)
 700             {
 701                 printf("FIT: Integral of fitted function: %g\n", integral);
 702                 if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn))
 703                 {
 704                     printf("FIT: Note that the constant term is not taken into account when "
 705                            "computing integral.\n");
 706                 }
 707             }
 708             /* Generate debug output */
 709             if (nullptr != fn_fitted)
 710             {
 711                 fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)", "Data (t)", oenv);
 712                 for (i = 0; (i < effnNparams(eFitFn)); i++)
 713                 {
 714                     fprintf(fp, "# fitparms[%d] = %g\n", i, fitparms[i]);
 715                 }
 716                 for (j = 0; (j < nfitpnts); j++)
 717                 {
 718                     real ttt = x0 ? x0[j] : dt * j;
 719                     fprintf(fp, "%10.5e  %10.5e  %10.5e\n", x[j], y[j], (lmcurves[eFitFn])(ttt, fitparms));
 720                 }
 721                 xvgrclose(fp);
 722             }
 723         }
 724     }
 725
 726     sfree(x);
 727     sfree(y);
 728     sfree(dy);
 729
 730     return integral;
 731 }
 732
 733 real fit_acf(int                     ncorr,
 734              int                     fitfn,
 735              const gmx_output_env_t* oenv,
 736              gmx_bool                bVerbose,
 737              real                    tbeginfit,
 738              real                    tendfit,
 739              real                    dt,
 740              real                    c1[],
 741              real*                   fit)
 742 {
 743     double   fitparm[3];
 744     double   tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
 745     real*    sig;
 746     int      i, j, jmax, nf_int;
 747     gmx_bool bPrint;
 748
 749     GMX_ASSERT(effnNparams(fitfn) < static_cast<int>(sizeof(fitparm) / sizeof(double)),
 750                "Fitting function with more than 3 parameters not supported!");
 751     bPrint = bVerbose || bDebugMode();
 752
 753     if (bPrint)
 754     {
 755         printf("COR:\n");
 756     }
 757
 758     if (tendfit <= 0)
 759     {
 760         tendfit = ncorr * dt;
 761     }
 762     nf_int = std::min(ncorr, static_cast<int>(tendfit / dt));
 763     sum    = print_and_integrate(debug, nf_int, dt, c1, nullptr, 1);
 764
 765     if (bPrint)
 766     {
 767         printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n", 0.0, dt * nf_int, sum);
 768         printf("COR: Relaxation times are computed as fit to an exponential:\n");
 769         printf("COR:   %s\n", effnDescription(fitfn));
 770         printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n",
 771                tbeginfit,
 772                std::min(ncorr * dt, tendfit));
 773     }
 774
 775     tStart = 0;
 776     if (bPrint)
 777     {
 778         printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
 779                "Fit from",
 780                "Integral",
 781                "Tail Value",
 782                "Sum (ps)",
 783                " a1 (ps)",
 784                (effnNparams(fitfn) >= 2) ? " a2 ()" : "",
 785                (effnNparams(fitfn) >= 3) ? " a3 (ps)" : "");
 786     }
 787
 788     snew(sig, ncorr);
 789
 790     if (tbeginfit > 0)
 791     {
 792         jmax = 3;
 793     }
 794     else
 795     {
 796         jmax = 1;
 797     }
 798     for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr * dt)); j++)
 799     {
 800         /* Estimate the correlation time for better fitting */
 801         c_start     = -1;
 802         ct_estimate = 0;
 803         for (i = 0; (i < ncorr) && (dt * i < tStart || c1[i] > 0); i++)
 804         {
 805             if (c_start < 0)
 806             {
 807                 if (dt * i >= tStart)
 808                 {
 809                     c_start     = c1[i];
 810                     ct_estimate = 0.5 * c1[i];
 811                 }
 812             }
 813             else
 814             {
 815                 ct_estimate += c1[i];
 816             }
 817         }
 818         if (c_start > 0)
 819         {
 820             ct_estimate *= dt / c_start;
 821         }
 822         else
 823         {
 824             /* The data is strange, so we need to choose somehting */
 825             ct_estimate = tendfit;
 826         }
 827         if (debug)
 828         {
 829             fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
 830         }
 831
 832         if (fitfn == effnEXPEXP)
 833         {
 834             fitparm[0] = 0.002 * ncorr * dt;
 835             fitparm[1] = 0.95;
 836             fitparm[2] = 0.2 * ncorr * dt;
 837         }
 838         else
 839         {
 840             /* Good initial guess, this increases the probability of convergence */
 841             fitparm[0] = ct_estimate;
 842             fitparm[1] = 1.0;
 843             fitparm[2] = 1.0;
 844         }
 845
 846         /* Generate more or less appropriate sigma's */
 847         for (i = 0; i < ncorr; i++)
 848         {
 849             sig[i] = sqrt(ct_estimate + dt * i);
 850         }
 851
 852         nf_int    = std::min(ncorr, static_cast<int>((tStart + 1e-4) / dt));
 853         sum       = print_and_integrate(debug, nf_int, dt, c1, nullptr, 1);
 854         tail_corr = do_lmfit(
 855                 ncorr, c1, sig, dt, nullptr, tStart, tendfit, oenv, bDebugMode(), fitfn, fitparm, 0, nullptr);
 856         sumtot = sum + tail_corr;
 857         if (fit && ((jmax == 1) || (j == 1)))
 858         {
 859             double mfp[3];
 860             for (i = 0; (i < 3); i++)
 861             {
 862                 mfp[i] = fitparm[i];
 863             }
 864             for (i = 0; (i < ncorr); i++)
 865             {
 866                 fit[i] = lmcurves[fitfn](i * dt, mfp);
 867             }
 868         }
 869         if (bPrint)
 870         {
 871             printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
 872             for (i = 0; (i < effnNparams(fitfn)); i++)
 873             {
 874                 printf(" %11.4e", fitparm[i]);
 875             }
 876             printf("\n");
 877         }
 878         tStart += tbeginfit;
 879     }
 880     sfree(sig);
 881
 882     return sumtot;
 883 }