2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
42 #include "gromacs/fileio/xvgr.h"
43 #include "gromacs/gmxana/gstat.h"
44 #include "gromacs/legacyheaders/typedefs.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/futil.h"
48 #include "gromacs/utility/smalloc.h"
50 const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3};
52 const char *s_ffn[effnNR+2] = {
53 NULL, "none", "exp", "aexp", "exp_exp", "vac",
54 "exp5", "exp7", "exp9", "erffit", NULL, NULL
56 /* We don't allow errest as a choice on the command line */
58 const char *longs_ffn[effnNR] = {
62 "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
63 "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
64 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5",
65 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
66 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
67 "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
68 "y = a2*ee(a1,x) + (1-a2)*ee(a2,x)"
71 extern gmx_bool mrqmin_new(real x[], real y[], real sig[], int ndata, real a[],
72 int ia[], int ma, real **covar, real **alpha, real *chisq,
73 void (*funcs)(real, real [], real *, real []),
76 static real myexp(real x, real A, real tau)
78 if ((A == 0) || (tau == 0))
85 void erffit (real x, real a[], real *y, real dyda[])
88 * y=(a1+a2)/2-(a1-a2)/2*erf((x-a3)/a4^2)
96 erfarg = (x-a[3])/(a[4]*a[4]);
97 erfarg2 = erfarg*erfarg;
98 erfval = gmx_erf(erfarg)/2;
99 derf = M_2_SQRTPI*(a[1]-a[2])/2*exp(-erfarg2)/(a[4]*a[4]);
100 *y = (a[1]+a[2])/2-(a[1]-a[2])*erfval;
101 dyda[1] = 1/2-erfval;
102 dyda[2] = 1/2+erfval;
104 dyda[4] = 2*derf*erfarg;
109 static void exp_one_parm(real x, real a[], real *y, real dyda[])
121 dyda[1] = x*e1/sqr(a[1]);
124 static void exp_two_parm(real x, real a[], real *y, real dyda[])
136 dyda[1] = x*a[2]*e1/sqr(a[1]);
140 static void exp_3_parm(real x, real a[], real *y, real dyda[])
144 * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
152 *y = a[2]*e1 + (1-a[2])*e2;
153 dyda[1] = x*a[2]*e1/sqr(a[1]);
155 dyda[3] = x*(1-a[2])*e2/sqr(a[3]);
158 static void exp_5_parm(real x, real a[], real *y, real dyda[])
162 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5
170 *y = a[1]*e1 + a[3]*e2 + a[5];
174 fprintf(debug, "exp_5_parm called: x = %10.3f y = %10.3f\n"
175 "a = ( %8.3f %8.3f %8.3f %8.3f %8.3f)\n",
176 x, *y, a[1], a[2], a[3], a[4], a[5]);
179 dyda[2] = x*e1/sqr(a[2]);
181 dyda[4] = x*e2/sqr(a[4]);
185 static void exp_7_parm(real x, real a[], real *y, real dyda[])
189 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
198 *y = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7];
201 dyda[2] = x*e1/sqr(a[2]);
203 dyda[4] = x*e2/sqr(a[4]);
205 dyda[6] = x*e3/sqr(a[6]);
209 static void exp_9_parm(real x, real a[], real *y, real dyda[])
213 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
223 *y = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7]*e4 + a[9];
226 dyda[2] = x*e1/sqr(a[2]);
228 dyda[4] = x*e2/sqr(a[4]);
230 dyda[6] = x*e3/sqr(a[6]);
232 dyda[8] = x*e4/sqr(a[8]);
236 static void vac_2_parm(real x, real a[], real *y, real dyda[])
240 * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
242 * = exp(-v) (cosh(wv) + 1/w sinh(wv))
247 * For tranverse current autocorrelation functions:
249 * a2 = 4 tau (eta/rho) k^2
253 double v, det, omega, omega2, em, ec, es;
261 omega = sqrt(omega2);
264 ec = em*0.5*(exp(omega*v)+exp(-omega*v));
265 es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
269 ec = em*cos(omega*v);
270 es = em*sin(omega*v)/omega;
273 dyda[2] = (v/det*ec+(v-1/det)*es)/(-2.0);
274 dyda[1] = (1-det)*v/a[1]*es;
279 dyda[2] = -v*v*em*(0.5+v/6);
280 dyda[1] = v*v/a[1]*em;
284 static void errest_3_parm(real x, real a[], real *y, real dyda[])
290 e1 = exp(-x/a[1]) - 1;
298 e2 = exp(-x/a[3]) - 1;
307 v1 = 2*a[1]*(e1*a[1]/x + 1);
308 v2 = 2*a[3]*(e2*a[3]/x + 1);
309 *y = a[2]*v1 + (1-a[2])*v2;
310 dyda[1] = 2* a[2] *(v1/a[1] + e1);
311 dyda[3] = 2*(1 - a[2])*(v2/a[3] + e2);
323 typedef void (*myfitfn)(real x, real a[], real *y, real dyda[]);
324 myfitfn mfitfn[effnNR] = {
325 exp_one_parm, exp_one_parm, exp_two_parm, exp_3_parm, vac_2_parm,
326 exp_5_parm, exp_7_parm, exp_9_parm, erffit, errest_3_parm
329 real fit_function(int eFitFn, real *parm, real x)
331 static real y, dum[8];
333 mfitfn[eFitFn](x, parm-1, &y, dum);
338 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
339 static gmx_bool lmfit_exp(int nfit, real x[], real y[], real dy[], real ftol,
340 real parm[], real dparm[], gmx_bool bVerbose,
343 real chisq, ochisq, alamda;
344 real *a, **covar, **alpha, *dum;
346 int i, j, ma, mfit, *lista, *ia;
348 if ((eFitFn < 0) || (eFitFn >= effnNR))
350 gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)",
351 effnNR-1, eFitFn, __FILE__, __LINE__);
354 ma = mfit = nfp_ffn[eFitFn]; /* number of parameters to fit */
361 for (i = 1; (i < ma+1); i++)
364 ia[i] = 1; /* fixed bug B.S.S 19/11 */
365 snew(covar[i], ma+1);
366 snew(alpha[i], ma+1);
372 fprintf(stderr, "Will keep parameters fixed during fit procedure: %d\n",
375 for (i = 0; i < ma; i++)
385 fprintf(debug, "%d parameter fit\n", mfit);
389 alamda = -1; /* Starting value */
391 for (i = 0; (i < mfit); i++)
399 fprintf(stderr, "%4s %10s %10s %10s %10s %10s %10s\n",
400 "Step", "chi^2", "Lambda", "A1", "A2", "A3", "A4");
405 /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
406 * &chisq,expfn[mfit-1],&alamda)
408 if (!mrqmin_new(x-1, y-1, dy-1, nfit, a, ia, ma, covar, alpha, &chisq,
409 mfitfn[eFitFn], &alamda))
416 fprintf(stderr, "%4d %10.5e %10.5e %10.5e",
417 j, chisq, alamda, a[1]);
420 fprintf(stderr, " %10.5e", a[2]);
424 fprintf(stderr, " %10.5e", a[3]);
428 fprintf(stderr, " %10.5e", a[4]);
430 fprintf(stderr, "\n");
433 bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
434 ((ochisq == chisq)));
436 while (bCont && (alamda != 0.0) && (j < 50));
439 fprintf(stderr, "\n");
442 /* Now get the covariance matrix out */
445 /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
446 * &chisq,expfn[mfit-1],&alamda)
448 if (!mrqmin_new(x-1, y-1, dy-1, nfit, a, ia, ma, covar, alpha, &chisq,
449 mfitfn[eFitFn], &alamda))
454 for (j = 0; (j < mfit); j++)
457 dparm[j] = covar[j+1][j+1];
460 for (i = 0; (i < ma+1); i++)
474 real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
475 real begintimefit, real endtimefit, const output_env_t oenv,
476 gmx_bool bVerbose, int eFitFn, real fitparms[], int fix)
481 int i, j, nparm, nfitpnts;
487 nparm = nfp_ffn[eFitFn];
490 fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm);
491 fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
492 eFitFn, begintimefit, endtimefit, dt);
500 for (i = 0; i < ndata; i++)
502 ttt = x0 ? x0[i] : dt*i;
503 if (ttt >= begintimefit && ttt <= endtimefit)
508 /* mrqmin does not like sig to be zero */
519 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n",
520 j, i, x[j], y[j], dy[j]);
527 if (nfitpnts < nparm)
529 fprintf(stderr, "Not enough data points for fitting!\n");
537 for (i = 0; (i < nparm); i++)
539 parm[i] = fitparms[i];
543 if (!lmfit_exp(nfitpnts, x, y, dy, ftol, parm, dparm, bVerbose, eFitFn, fix))
545 fprintf(stderr, "Fit failed!\n");
549 /* Compute the integral from begintimefit */
552 integral = (parm[0]*myexp(begintimefit, parm[1], parm[0]) +
553 parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));
557 integral = parm[0]*myexp(begintimefit, parm[1], parm[0]);
561 integral = parm[0]*myexp(begintimefit, 1, parm[0]);
565 gmx_fatal(FARGS, "nparm = %d in file %s, line %d",
566 nparm, __FILE__, __LINE__);
569 /* Generate THE output */
572 fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts);
573 fprintf(stderr, "FIT: %21s%21s%21s\n",
574 "parm0 ", "parm1 (ps) ", "parm2 (ps) ");
575 fprintf(stderr, "FIT: ------------------------------------------------------------\n");
576 fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
577 parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);
578 fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
579 begintimefit, integral);
581 sprintf(buf, "test%d.xvg", nfitpnts);
582 fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);
583 fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n",
584 parm[0], parm[1], parm[2]);
585 for (j = 0; (j < nfitpnts); j++)
587 ttt = x0 ? x0[j] : dt*j;
588 fprintf(fp, "%10.5e %10.5e %10.5e\n",
589 ttt, c1[j], fit_function(eFitFn, parm, ttt));
596 for (i = 0; (i < nparm); i++)
598 fitparms[i] = parm[i];
612 void do_expfit(int ndata, real c1[], real dt, real begintimefit, real endtimefit)
616 real aa, bb, saa, sbb, A, tau, dA, dtau;
618 fprintf(stderr, "Will fit data from %g (ps) to %g (ps).\n",
619 begintimefit, endtimefit);
621 snew(x, ndata); /* allocate the maximum necessary space */
626 for (i = 0; (i < ndata); i++)
628 if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) )
633 fprintf(stderr, "n= %d, i= %d, x= %g, y= %g\n", n, i, x[n], y[n]);
637 fprintf(stderr, "# of data points used in the fit is : %d\n\n", n);
638 expfit(n, x, y, Dy, &aa, &saa, &bb, &sbb);
644 fprintf(stderr, "Fitted to y=exp(a+bx):\n");
645 fprintf(stderr, "a = %10.5f\t b = %10.5f", aa, bb);
646 fprintf(stderr, "\n");
647 fprintf(stderr, "Fitted to y=Aexp(-x/tau):\n");
648 fprintf(stderr, "A = %10.5f\t tau = %10.5f\n", A, tau);
649 fprintf(stderr, "dA = %10.5f\t dtau = %10.5f\n", dA, dtau);
653 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
656 real *w, *ly, A, SA, B, SB;
658 real sum, xbar, ybar, Sxx, Sxy, wr2, chi2, gamma, Db;
665 #define sqr(x) ((x)*(x))
671 /* Calculate weights and values of ln(y). */
672 for (i = 0; (i < n); i++)
674 w[i] = sqr(y[i]/Dy[i]);
678 /* The weighted averages of x and y: xbar and ybar. */
682 for (i = 0; (i < n); i++)
691 /* The centered product sums Sxx and Sxy, and hence A and B. */
694 for (i = 0; (i < n); i++)
696 Sxx += w[i]*sqr(x[i]-xbar);
697 Sxy += w[i]*(x[i]-xbar)*(ly[i]-ybar);
702 /* Chi-squared (chi2) and gamma. */
705 for (i = 0; (i < n); i++)
707 wr2 = w[i]*sqr(ly[i]-A-B*x[i]);
709 gamma += wr2*(x[i]-xbar);
712 /* Refined values of A and B. Also SA and SB. */
713 Db = -ONEP5*gamma/Sxx;
715 A -= ONEP5*chi2/sum-xbar*Db;
716 SB = sqrt((chi2/(n-2))/Sxx);
717 SA = SB*sqrt(Sxx/sum+sqr(xbar));