src/gromacs/correlationfunctions/expfit.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  * Copyright (c) 2013,2014, 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 <math.h>
  50 #include <string.h>
  51
  52 #include "external/lmfit/lmcurve.h"
  53
  54 #include "gromacs/correlationfunctions/integrate.h"
  55 #include "gromacs/fileio/xvgr.h"
  56 #include "gromacs/legacyheaders/macros.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/futil.h"
  60 #include "gromacs/utility/real.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 /*! \brief Number of parameters for each fitting function */
  64 static const int   nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3, 6};
  65
  66 const char        *s_ffn[effnNR+2] = {
  67     NULL, "none", "exp", "aexp", "exp_exp", "vac",
  68     "exp5", "exp7", "exp9", "erffit", "effnPRES", NULL, NULL
  69 };
  70 /* We don't allow errest as a choice on the command line */
  71
  72 /*! \brief Long description for each fitting function type */
  73 static const char *longs_ffn[effnNR] = {
  74     "no fit",
  75     "y = exp(-x/a1)",
  76     "y = a2 exp(-x/a1)",
  77     "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
  78     "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
  79     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5",
  80     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
  81     "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
  82     "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
  83     "y = a2 *2*a1*((exp(-x/a1)-1)*(a1/x)+1)+(1-a2)*2*a3*((exp(-x/a3)-1)*(a3/x)+1)",
  84     "y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6)"
  85 };
  86
  87 int effnNparams(int effn)
  88 {
  89     if ((0 <= effn) && (effn < effnNR))
  90     {
  91         return nfp_ffn[effn];
  92     }
  93     else
  94     {
  95         return -1;
  96     }
  97 }
  98
  99 const char *effnDescription(int effn)
 100 {
 101     if ((0 <= effn) && (effn < effnNR))
 102     {
 103         return longs_ffn[effn];
 104     }
 105     else
 106     {
 107         return NULL;
 108     }
 109 }
 110
 111 int sffn2effn(const char **sffn)
 112 {
 113     int eFitFn, i;
 114
 115     eFitFn = 0;
 116     for (i = 0; i < effnNR; i++)
 117     {
 118         if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
 119         {
 120             eFitFn = i;
 121         }
 122     }
 123
 124     return eFitFn;
 125 }
 126
 127 /*! \brief Compute exponential function A exp(-x/tau) */
 128 static double myexp(double x, double A, double tau)
 129 {
 130     if ((A == 0) || (tau == 0))
 131     {
 132         return 0;
 133     }
 134     return A*exp(-x/tau);
 135 }
 136
 137 /*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
 138 static double lmc_erffit (double x, const double *a)
 139 {
 140     double erfarg;
 141
 142     erfarg  = (x-a[2])/(a[3]*a[3]);
 143     return (a[0]+a[1])/2-(a[0]-a[1])*gmx_erf(erfarg);
 144 }
 145
 146 /*! \brief Compute y = exp(-x/a0) */
 147 static double lmc_exp_one_parm(double x, const double *a)
 148 {
 149     return exp(-x/a[0]);
 150 }
 151
 152 /*! \brief Compute y = a1 exp(-x/a0) */
 153 static double lmc_exp_two_parm(double x, const double *a)
 154 {
 155     return a[1]*exp(-x/a[0]);
 156 }
 157
 158 /*! \brief Compute y = a1 exp(-x/a0) + (1-a1) exp(-x/a2) */
 159 static double lmc_exp_3_parm(double x, const double *a)
 160 {
 161     double e1, e2;
 162
 163     e1      = exp(-x/a[0]);
 164     e2      = exp(-x/a[2]);
 165     return a[1]*e1 + (1-a[1])*e2;
 166 }
 167
 168 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 */
 169 static double lmc_exp_5_parm(double x, const double *a)
 170 {
 171     double e1, e2;
 172
 173     e1      = exp(-x/a[1]);
 174     e2      = exp(-x/a[3]);
 175     return a[0]*e1 + a[2]*e2 + a[4];
 176 }
 177
 178 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 */
 179 static double lmc_exp_7_parm(double x, const double *a)
 180 {
 181     double e1, e2, e3;
 182
 183     e1      = exp(-x/a[1]);
 184     e2      = exp(-x/a[3]);
 185     e3      = exp(-x/a[5]);
 186     return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6];
 187 }
 188
 189 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 exp(-x/a7) + a8 */
 190 static double lmc_exp_9_parm(double x, const double *a)
 191 {
 192     double e1, e2, e3, e4;
 193
 194     e1      = exp(-x/a[1]);
 195     e2      = exp(-x/a[3]);
 196     e3      = exp(-x/a[5]);
 197     e4      = exp(-x/a[7]);
 198     return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8];
 199 }
 200
 201 /*! \brief Compute y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6) */
 202 static double effnPRES_fun(double x, const double *a)
 203 {
 204     double term1, term2, term3;
 205
 206     if (a[0] != 0)
 207     {
 208         term3 = a[0] * exp(-pow((x/a[4]), a[5]));
 209     }
 210     else
 211     {
 212         term3 = 0;
 213     }
 214
 215     term1 = 1-a[0];
 216     if (term1 != 0)
 217     {
 218         term2 = exp(-pow((x/a[2]), a[3])) * cos(x*a[1]);
 219     }
 220     else
 221     {
 222         term2 = 0;
 223     }
 224
 225     return term1*term2 + term3;
 226 }
 227
 228 /*! \brief Compute vac function */
 229 static double lmc_vac_2_parm(double x, const double *a)
 230 {
 231     /* Fit to function
 232      *
 233      * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 234      *
 235      *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 236      *
 237      *    v = x/(2 a1)
 238      *    w = sqrt(1 - a2)
 239      *
 240      *    For tranverse current autocorrelation functions:
 241      *       a1 = tau
 242      *       a2 = 4 tau (eta/rho) k^2
 243      *
 244      */
 245
 246     double y, v, det, omega, omega2, em, ec, es;
 247
 248     v   = x/(2*a[0]);
 249     det = 1 - a[1];
 250     em  = exp(-v);
 251     if (det != 0)
 252     {
 253         omega2  = fabs(det);
 254         omega   = sqrt(omega2);
 255         if (det > 0)
 256         {
 257             ec = em*0.5*(exp(omega*v)+exp(-omega*v));
 258             es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
 259         }
 260         else
 261         {
 262             ec = em*cos(omega*v);
 263             es = em*sin(omega*v)/omega;
 264         }
 265         y      = ec + es;
 266     }
 267     else
 268     {
 269         y      = (1+v)*em;
 270     }
 271     return y;
 272 }
 273
 274 /*! \brief Compute error estimate */
 275 static double lmc_errest_3_parm(double x, const double *a)
 276 {
 277     /*
 278        e1 = (exp(-x/a1) - 1)
 279        e2 = (exp(-x/a3) - 1)
 280        v1=  2*a1 * (e1*a1/x + 1)
 281        v2 = 2*a3 * (e2*a3/x + 1)
 282        fun = a2*v1 + (1 - a2) * v2
 283      */
 284     double e1, e2, v1, v2;
 285
 286     if (a[0] != 0)
 287     {
 288         e1 = exp(-x/a[0]) - 1;
 289     }
 290     else
 291     {
 292         e1 = 0;
 293     }
 294     if (a[2] != 0)
 295     {
 296         e2 = exp(-x/a[2]) - 1;
 297     }
 298     else
 299     {
 300         e2 = 0;
 301     }
 302
 303     if (x > 0)
 304     {
 305         v1      = 2*a[0]*(e1*a[0]/x + 1);
 306         v2      = 2*a[2]*(e2*a[2]/x + 1);
 307         return a[1]*v1 + (1-a[1])*v2;
 308     }
 309     else
 310     {
 311         return 0;
 312     }
 313 }
 314
 315 /*! \brief function type for passing to fitting routine */
 316 typedef double (*t_lmcurve)(double x, const double *a);
 317
 318 /*! \brief array of fitting functions corresponding to the pre-defined types */
 319 t_lmcurve lmcurves[effnNR] = {
 320     lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
 321     lmc_exp_3_parm, lmc_vac_2_parm,
 322     lmc_exp_5_parm,   lmc_exp_7_parm,
 323     lmc_exp_9_parm, lmc_erffit,  lmc_errest_3_parm, effnPRES_fun
 324 };
 325
 326 double fit_function(int eFitFn, double *parm, double x)
 327 {
 328     // TODO: check that eFitFn is within bounds
 329     return lmcurves[eFitFn](x, parm);
 330 }
 331
 332 /*! \brief lmfit_exp supports up to 3 parameter fitting of exponential functions */
 333 static gmx_bool lmfit_exp(int nfit, double x[], double y[],
 334                           real ftol, int maxiter,
 335                           double parm[], gmx_bool bVerbose,
 336                           int eFitFn)
 337 {
 338     double             chisq, ochisq;
 339     gmx_bool           bCont;
 340     int                i, j;
 341     lm_control_struct *control;
 342     lm_status_struct  *status;
 343
 344     if ((eFitFn < 0) || (eFitFn >= effnNR))
 345     {
 346         gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)",
 347                   effnNR-1, eFitFn, __FILE__, __LINE__);
 348     }
 349
 350     snew(control, 1);
 351     control->ftol       = ftol;
 352     control->xtol       = ftol;
 353     control->gtol       = ftol;
 354     control->epsilon    = 0.1;
 355     control->stepbound  = 100;
 356     control->patience   = maxiter;
 357     control->scale_diag = 1;
 358     control->msgfile    = stdout;
 359     control->verbosity  = (bVerbose ? 1 : 0);
 360     control->n_maxpri   = nfp_ffn[eFitFn];
 361     control->m_maxpri   = 0;
 362     snew(status, 1);
 363     /* Initial params */
 364     chisq  = 1e12;
 365     j      = 0;
 366     if (bVerbose)
 367     {
 368         fprintf(stderr, "%4s  %10s  Parameters\n",
 369                 "Step", "chi^2");
 370     }
 371     do
 372     {
 373         ochisq = chisq;
 374
 375         lmcurve(nfp_ffn[eFitFn], parm, nfit, x, y,
 376                 lmcurves[eFitFn], control, status);
 377         chisq = sqr(status->fnorm);
 378         if (bVerbose)
 379         {
 380             int mmm;
 381             fprintf(stderr, "%4d  %10.5e", j, chisq);
 382             for (mmm = 0; (mmm < nfp_ffn[eFitFn]); mmm++)
 383             {
 384                 fprintf(stderr, "  %10.5e", parm[mmm]);
 385             }
 386             fprintf(stderr, "\n");
 387         }
 388         j++;
 389         bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
 390                  ((ochisq == chisq)));
 391     }
 392     while (bCont && (j < maxiter));
 393     if (bVerbose)
 394     {
 395         fprintf(stderr, "\n");
 396     }
 397
 398     sfree(control);
 399     sfree(status);
 400
 401     return TRUE;
 402 }
 403
 404 /*! \brief See description in header file. */
 405 real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
 406               real begintimefit, real endtimefit, const output_env_t oenv,
 407               gmx_bool bVerbose, int eFitFn, double fitparms[], int fix)
 408 {
 409     FILE    *fp;
 410     char     buf[32];
 411
 412     int      i, j, nparm, nfitpnts;
 413     double   integral, ttt;
 414     double  *parm, *dparm;
 415     double  *x, *y, *dy;
 416 #ifdef GMX_DOUBLE
 417     real     ftol    = 1e-16;
 418 #else
 419     real     ftol    = 1e-8;
 420 #endif
 421     int      maxiter = 50;
 422
 423     if (0 != fix)
 424     {
 425         fprintf(stderr, "Using fixed parameters in curve fitting is not working anymore\n");
 426     }
 427     nparm = nfp_ffn[eFitFn];
 428     if (debug)
 429     {
 430         fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm);
 431         fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
 432                 eFitFn, begintimefit, endtimefit, dt);
 433     }
 434
 435     snew(x, ndata);
 436     snew(y, ndata);
 437     snew(dy, ndata);
 438
 439     j = 0;
 440     for (i = 0; i < ndata; i++)
 441     {
 442         ttt = x0 ? x0[i] : dt*i;
 443         if (ttt >= begintimefit && ttt <= endtimefit)
 444         {
 445             x[j] = ttt;
 446             y[j] = c1[i];
 447
 448             /* mrqmin does not like sig to be zero */
 449             if (sig[i] < 1.0e-7)
 450             {
 451                 dy[j] = 1.0e-7;
 452             }
 453             else
 454             {
 455                 dy[j] = sig[i];
 456             }
 457             if (debug)
 458             {
 459                 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n",
 460                         j, i, x[j], y[j], dy[j]);
 461             }
 462             j++;
 463         }
 464     }
 465     nfitpnts = j;
 466     integral = 0;
 467     if (nfitpnts < nparm)
 468     {
 469         fprintf(stderr, "Not enough data points for fitting!\n");
 470     }
 471     else
 472     {
 473         snew(parm, nparm);
 474         snew(dparm, nparm);
 475         if (fitparms)
 476         {
 477             for (i = 0; (i < nparm); i++)
 478             {
 479                 parm[i] = fitparms[i];
 480             }
 481         }
 482
 483         if (!lmfit_exp(nfitpnts, x, y, ftol, maxiter, parm, bVerbose, eFitFn))
 484         {
 485             fprintf(stderr, "Fit failed!\n");
 486         }
 487         else if (nparm <= 3)
 488         {
 489             /* Compute the integral from begintimefit */
 490             if (nparm == 3)
 491             {
 492                 integral = (parm[0]*myexp(begintimefit, parm[1],  parm[0]) +
 493                             parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));
 494             }
 495             else if (nparm == 2)
 496             {
 497                 integral = parm[0]*myexp(begintimefit, parm[1],  parm[0]);
 498             }
 499             else if (nparm == 1)
 500             {
 501                 integral = parm[0]*myexp(begintimefit, 1,  parm[0]);
 502             }
 503             else
 504             {
 505                 gmx_fatal(FARGS, "nparm = %d in file %s, line %d",
 506                           nparm, __FILE__, __LINE__);
 507             }
 508
 509             /* Generate THE output */
 510             if (bVerbose)
 511             {
 512                 fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts);
 513                 fprintf(stderr, "FIT: %21s%21s%21s\n",
 514                         "parm0     ", "parm1 (ps)   ", "parm2 (ps)    ");
 515                 fprintf(stderr, "FIT: ------------------------------------------------------------\n");
 516                 fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
 517                         parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);
 518                 fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
 519                         begintimefit, integral);
 520
 521                 sprintf(buf, "test%d.xvg", nfitpnts);
 522                 fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);
 523                 fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n",
 524                         parm[0], parm[1], parm[2]);
 525                 for (j = 0; (j < nfitpnts); j++)
 526                 {
 527                     ttt = x0 ? x0[j] : dt*j;
 528                     fprintf(fp, "%10.5e  %10.5e  %10.5e\n",
 529                             ttt, c1[j],
 530                             lmcurves[eFitFn](ttt, parm));
 531                 }
 532                 xvgrclose(fp);
 533             }
 534         }
 535         if (fitparms)
 536         {
 537             for (i = 0; (i < nparm); i++)
 538             {
 539                 fitparms[i] = parm[i];
 540             }
 541         }
 542         sfree(parm);
 543         sfree(dparm);
 544     }
 545
 546     sfree(x);
 547     sfree(y);
 548     sfree(dy);
 549
 550     return integral;
 551 }
 552
 553 /*! See description in header file. */
 554 real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
 555              real tbeginfit, real tendfit, real dt, real c1[], real *fit)
 556 {
 557     double      fitparm[3];
 558     double      tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
 559     real       *sig;
 560     int         i, j, jmax, nf_int;
 561     gmx_bool    bPrint;
 562
 563     bPrint = bVerbose || bDebugMode();
 564
 565     if (bPrint)
 566     {
 567         printf("COR:\n");
 568     }
 569
 570     if (tendfit <= 0)
 571     {
 572         tendfit = ncorr*dt;
 573     }
 574     nf_int = min(ncorr, (int)(tendfit/dt));
 575     sum    = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
 576
 577     /* Estimate the correlation time for better fitting */
 578     ct_estimate = 0.5*c1[0];
 579     for (i = 1; (i < ncorr) && (c1[i] > 0); i++)
 580     {
 581         ct_estimate += c1[i];
 582     }
 583     ct_estimate *= dt/c1[0];
 584
 585     if (bPrint)
 586     {
 587         printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
 588                0.0, dt*nf_int, sum);
 589         printf("COR: Relaxation times are computed as fit to an exponential:\n");
 590         printf("COR:   %s\n", longs_ffn[fitfn]);
 591         printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, min(ncorr*dt, tendfit));
 592     }
 593
 594     tStart = 0;
 595     if (bPrint)
 596     {
 597         printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
 598                "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
 599                (nfp_ffn[fitfn] >= 2) ? " a2 ()" : "",
 600                (nfp_ffn[fitfn] >= 3) ? " a3 (ps)" : "");
 601     }
 602
 603     snew(sig, ncorr);
 604
 605     if (tbeginfit > 0)
 606     {
 607         jmax = 3;
 608     }
 609     else
 610     {
 611         jmax = 1;
 612     }
 613     for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++)
 614     {
 615         /* Estimate the correlation time for better fitting */
 616         c_start     = -1;
 617         ct_estimate = 0;
 618         for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++)
 619         {
 620             if (c_start < 0)
 621             {
 622                 if (dt*i >= tStart)
 623                 {
 624                     c_start     = c1[i];
 625                     ct_estimate = 0.5*c1[i];
 626                 }
 627             }
 628             else
 629             {
 630                 ct_estimate += c1[i];
 631             }
 632         }
 633         if (c_start > 0)
 634         {
 635             ct_estimate *= dt/c_start;
 636         }
 637         else
 638         {
 639             /* The data is strange, so we need to choose somehting */
 640             ct_estimate = tendfit;
 641         }
 642         if (debug)
 643         {
 644             fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
 645         }
 646
 647         if (fitfn == effnEXP3)
 648         {
 649             fitparm[0] = 0.002*ncorr*dt;
 650             fitparm[1] = 0.95;
 651             fitparm[2] = 0.2*ncorr*dt;
 652         }
 653         else
 654         {
 655             /* Good initial guess, this increases the probability of convergence */
 656             fitparm[0] = ct_estimate;
 657             fitparm[1] = 1.0;
 658             fitparm[2] = 1.0;
 659         }
 660
 661         /* Generate more or less appropriate sigma's */
 662         for (i = 0; i < ncorr; i++)
 663         {
 664             sig[i] = sqrt(ct_estimate+dt*i);
 665         }
 666
 667         nf_int    = min(ncorr, (int)((tStart+1e-4)/dt));
 668         sum       = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
 669         tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv,
 670                              bDebugMode(), fitfn, fitparm, 0);
 671         sumtot = sum+tail_corr;
 672         if (fit && ((jmax == 1) || (j == 1)))
 673         {
 674             double mfp[3];
 675             for (i = 0; (i < 3); i++)
 676             {
 677                 mfp[i] = fitparm[i];
 678             }
 679             for (i = 0; (i < ncorr); i++)
 680             {
 681                 fit[i] = lmcurves[fitfn](i*dt, mfp);
 682             }
 683         }
 684         if (bPrint)
 685         {
 686             printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
 687             for (i = 0; (i < nfp_ffn[fitfn]); i++)
 688             {
 689                 printf(" %11.4e", fitparm[i]);
 690             }
 691             printf("\n");
 692         }
 693         tStart += tbeginfit;
 694     }
 695     sfree(sig);
 696
 697     return sumtot;
 698 }