6 #include <gsl/gsl_rng.h>
7 #include <gsl/gsl_randist.h>
8 #include <gsl/gsl_vector.h>
9 #include <gsl/gsl_blas.h>
10 #include <gsl/gsl_multimin.h>
11 #include <gsl/gsl_multifit_nlin.h>
12 #include <gsl/gsl_sf.h>
13 #include <gsl/gsl_version.h>
28 /* The first few sections of this file contain functions that were adopted,
29 * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
30 * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
31 * This is also the case with the function eq10v2().
33 * The parts menetioned in the previous paragraph were contributed under a BSD license.
37 /* This first part is from complex.c which I recieved from Omer Markowitch.
40 * ------------- from complex.c ------------- */
42 /* Complexation of a paired number (r,i) */
43 static gem_complex
gem_cmplx(double r
, double i
)
51 /* Complexation of a real number, x */
52 static gem_complex
gem_c(double x
)
60 /* Real and Imaginary part of a complex number z -- Re (z) and Im (z) */
61 static double gem_Re(gem_complex z
) {return z
.r
;}
62 static double gem_Im(gem_complex z
) {return z
.i
;}
64 /* Magnitude of a complex number z */
65 static double gem_cx_abs(gem_complex z
) { return (sqrt(z
.r
*z
.r
+z
.i
*z
.i
)); }
67 /* Addition of two complex numbers z1 and z2 */
68 static gem_complex
gem_cxadd(gem_complex z1
, gem_complex z2
)
76 /* Addition of a complex number z1 and a real number r */
77 static gem_complex
gem_cxradd(gem_complex z
, double r
)
85 /* Subtraction of two complex numbers z1 and z2 */
86 static gem_complex
gem_cxsub(gem_complex z1
, gem_complex z2
)
94 /* Multiplication of two complex numbers z1 and z2 */
95 static gem_complex
gem_cxmul(gem_complex z1
, gem_complex z2
)
98 value
.r
= z1
.r
*z2
.r
-z1
.i
*z2
.i
;
99 value
.i
= z1
.r
*z2
.i
+z1
.i
*z2
.r
;
103 /* Square of a complex number z */
104 static gem_complex
gem_cxsq(gem_complex z
)
107 value
.r
= z
.r
*z
.r
-z
.i
*z
.i
;
108 value
.i
= z
.r
*z
.i
*2.;
112 /* multiplication of a complex number z and a real number r */
113 static gem_complex
gem_cxrmul(gem_complex z
, double r
)
121 /* Division of two complex numbers z1 and z2 */
122 static gem_complex
gem_cxdiv(gem_complex z1
, gem_complex z2
)
126 num
= z2
.r
*z2
.r
+z2
.i
*z2
.i
;
129 fprintf(stderr
, "ERROR in gem_cxdiv function\n");
131 value
.r
= (z1
.r
*z2
.r
+z1
.i
*z2
.i
)/num
; value
.i
= (z1
.i
*z2
.r
-z1
.r
*z2
.i
)/num
;
135 /* division of a complex z number by a real number x */
136 static gem_complex
gem_cxrdiv(gem_complex z
, double r
)
144 /* division of a real number r by a complex number x */
145 static gem_complex
gem_rxcdiv(double r
, gem_complex z
)
149 f
= r
/(z
.r
*z
.r
+z
.i
*z
.i
);
155 /* Integer power of a complex number z -- z^x */
156 static gem_complex
gem_cxintpow(gem_complex z
, int x
)
167 value
= gem_cxmul(value
, z
);
174 value
= gem_cxdiv(value
, z
);
183 /* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
184 static gem_complex
gem_cxdexp(gem_complex z
)
189 value
.r
= exp_z_r
*cos(z
.i
);
190 value
.i
= exp_z_r
*sin(z
.i
);
194 /* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z), */
195 /* where -PI < Arg(z) < PI */
196 static gem_complex
gem_cxlog(gem_complex z
)
200 mag2
= z
.r
*z
.r
+z
.i
*z
.i
;
202 fprintf(stderr
, "ERROR in gem_cxlog func\n");
204 value
.r
= log(sqrt(mag2
));
211 value
.i
= atan2(z
.i
, z
.r
);
216 /* Square root of a complex number z = |z| exp(I*the) -- z^(1/2) */
217 /* z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */
218 /* where 0 < the < 2*PI */
219 static gem_complex
gem_cxdsqrt(gem_complex z
)
224 value
.r
= sqrt(fabs((sq
+z
.r
)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */
225 value
.i
= sqrt(fabs((sq
-z
.r
)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */
232 /* square root of a real number r */
233 static gem_complex
gem_cxrsqrt(double r
) {
236 return(gem_cmplx(0, sqrt(-r
)));
240 return(gem_c(sqrt(r
)));
244 /* Complex power of a complex number z1^z2 */
245 static gem_complex
gem_cxdpow(gem_complex z1
, gem_complex z2
)
248 value
=gem_cxdexp(gem_cxmul(gem_cxlog(z1
),z2
));
252 /* Print out a complex number z as z: z.r, z.i */
253 static void gem_cxprintf(gem_complex z
) { fprintf(stdout
, "z: %lg + %lg_i\n", z
.r
, z
.i
); }
254 /* ------------ end of complex.c ------------ */
256 /* This next part was derived from cubic.c, also received from Omer Markovitch.
257 * ------------- from cubic.c ------------- */
259 /* Solver for a cubic equation: x^3-a*x^2+b*x-c=0 */
260 static void gem_solve(gem_complex
*al
, gem_complex
*be
, gem_complex
*gam
,
261 double a
, double b
, double c
)
263 double t1
, t2
, two_3
, temp
;
264 gem_complex ctemp
, ct3
;
266 two_3
=pow(2., 1./3.); t1
= -a
*a
+3.*b
; t2
= 2.*a
*a
*a
-9.*a
*b
+27.*c
;
267 temp
= 4.*t1
*t1
*t1
+t2
*t2
;
269 ctemp
= gem_cmplx(temp
,0.); ctemp
= gem_cxadd(gem_cmplx(t2
,0.),gem_cxdsqrt(ctemp
));
270 ct3
= gem_cxdpow(ctemp
, gem_cmplx(1./3.,0.));
272 ctemp
= gem_rxcdiv(-two_3
*t1
/3., ct3
);
273 ctemp
= gem_cxadd(ctemp
, gem_cxrdiv(ct3
, 3.*two_3
));
275 *gam
= gem_cxadd(gem_cmplx(a
/3.,0.), ctemp
);
277 ctemp
=gem_cxmul(gem_cxsq(*gam
), gem_cxsq(gem_cxsub(*gam
, gem_cmplx(a
,0.))));
278 ctemp
=gem_cxadd(ctemp
, gem_cxmul(gem_cmplx(-4.*c
,0.), *gam
));
279 ctemp
= gem_cxdiv(gem_cxdsqrt(ctemp
), *gam
);
280 *al
= gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a
,0.), *gam
),ctemp
),0.5);
281 *be
= gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a
,0.), *gam
),ctemp
),0.5);
284 /* ------------ end of cubic.c ------------ */
286 /* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
287 * ------------- from [cr]error.c ------------- */
289 /************************************************************/
290 /* Real valued error function and related functions */
291 /* x, y : real variables */
292 /* erf(x) : error function */
293 /* erfc(x) : complementary error function */
294 /* omega(x) : exp(x*x)*erfc(x) */
295 /* W(x,y) : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
296 /************************************************************/
298 /*---------------------------------------------------------------------------*/
299 /* Utilzed the series approximation of erf(x) */
300 /* Relative error=|err(x)|/erf(x)<EPS */
301 /* Handbook of Mathematical functions, Abramowitz, p 297 */
302 /* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
303 /*---------------------------------------------------------------------------*/
304 /* This gives erfc(z) correctly only upto >10-15 */
306 static double gem_erf(double x
)
308 double n
,sum
,temp
,exp2
,xm
,x2
,x4
,x6
,x8
,x10
,x12
;
321 for(n
=1.; n
<=2000.; n
+=1.){
322 temp
*= 2.*x2
/(2.*n
+1.);
324 if(fabs(temp
/sum
)<1.E
-16)
329 fprintf(stderr
, "In Erf calc - iteration exceeds %lg\n",n
);
334 /* from the asymptotic expansion of experfc(x) */
335 sum
= (1. - 0.5/x2
+ 0.75/x4
336 - 1.875/x6
+ 6.5625/x8
337 - 29.53125/x10
+ 162.421875/x12
)
339 sum
*=exp2
; /* now sum is erfc(x) */
342 return x
>=0.0 ? sum
: -sum
;
345 /* Result --> Alex's code for erfc and experfc behaves better than this */
346 /* complementray error function. Returns 1.-erf(x) */
347 static double gem_erfc(double x
)
353 ans
= t
* exp(-z
*z
- 1.26551223 +
362 t
*0.17087277)))))))));
364 return x
>=0.0 ? ans
: 2.0-ans
;
367 /* omega(x)=exp(x^2)erfc(x) */
368 static double gem_omega(double x
)
370 double xm
, ans
, xx
, x4
, x6
, x8
, x10
, x12
;
380 ans
= exp(xx
)*gem_erfc(x
);
382 /* Asymptotic expansion */
383 ans
= (1. - 0.5/xx
+ 0.75/x4
- 1.875/x6
+ 6.5625/x8
- 29.53125/x10
+ 162.421875/x12
) / sPI
/x
;
388 /* W(x,y)=exp(-x^2)*omega(x+y)=exp(2xy+y^2)*erfc(x+y) */
389 static double gem_W(double x
, double y
){ return(exp(-x
*x
)*gem_omega(x
+y
)); }
391 /**************************************************************/
392 /* Complex error function and related functions */
393 /* x, y : real variables */
394 /* z : complex variable */
395 /* cerf(z) : error function */
396 /* comega(z): exp(z*z)*cerfc(z) */
397 /* W(x,z) : exp(-x*x)*comega(x+z)=exp(2*x*z+z^2)*cerfc(x+z) */
398 /**************************************************************/
399 static gem_complex
gem_cerf(gem_complex z
)
403 double sumr
,sumi
,n
,n2
,f
,temp
,temp1
;
404 double x2
,cos_2xy
,sin_2xy
,cosh_2xy
,sinh_2xy
,cosh_ny
,sinh_ny
;
411 cos_2xy
= cos(2.*x
*y
);
412 sin_2xy
= sin(2.*x
*y
);
413 cosh_2xy
= cosh(2.*x
*y
);
414 sinh_2xy
= sinh(2.*x
*y
);
416 for(n
=1.0,temp
=0.; n
<=2000.; n
+=1.0)
421 f
= exp(-n2
/4.)/(n2
+4.*x2
);
422 sumr
+= (2.*x
- 2.*x
*cosh_ny
*cos_2xy
+n
*sinh_ny
*sin_2xy
)*f
;
423 sumi
+= (2.*x
*cosh_ny
*sin_2xy
+ n
*sinh_ny
*cos_2xy
)*f
;
424 temp1
= sqrt(sumr
*sumr
+sumi
*sumi
);
425 if(fabs((temp1
-temp
)/temp1
)<1.E
-16) {
432 fprintf(stderr
, "iteration exceeds %lg\n",n
);
443 value
.r
= gem_erf(x
) + (f
*(1.-cos_2xy
) + sumr
)*exp(-x2
);
444 value
.i
= (f
*sin_2xy
+sumi
)*exp(-x2
);
448 /*---------------------------------------------------------------------------*/
449 /* Utilzed the series approximation of erf(z=x+iy) */
450 /* Relative error=|err(z)|/|erf(z)|<EPS */
451 /* Handbook of Mathematical functions, Abramowitz, p 299 */
452 /* comega(z=x+iy)=cexp(z^2)*cerfc(z) */
453 /*---------------------------------------------------------------------------*/
454 static gem_complex
gem_comega(gem_complex z
)
458 double sumr
,sumi
,n
,n2
,f
,temp
,temp1
;
459 double x2
, y2
,cos_2xy
,sin_2xy
,cosh_2xy
,sinh_2xy
,cosh_ny
,sinh_ny
,exp_y2
;
467 cos_2xy
= cos(2.*x
*y
);
468 sin_2xy
= sin(2.*x
*y
);
469 cosh_2xy
= cosh(2.*x
*y
);
470 sinh_2xy
= sinh(2.*x
*y
);
473 for(n
=1.0, temp
=0.; n
<=2000.; n
+=1.0){
477 f
= exp(-n2
/4.)/(n2
+4.*x2
);
478 /* if(f<1.E-200) break; */
479 sumr
+= (2.*x
- 2.*x
*cosh_ny
*cos_2xy
+ n
*sinh_ny
*sin_2xy
)*f
;
480 sumi
+= (2.*x
*cosh_ny
*sin_2xy
+ n
*sinh_ny
*cos_2xy
)*f
;
481 temp1
= sqrt(sumr
*sumr
+sumi
*sumi
);
482 if(fabs((temp1
-temp
)/temp1
)<1.E
-16) {
488 fprintf(stderr
, "iteration exceeds %lg\n",n
);
498 value
.r
= gem_omega(x
)-(f
*(1.-cos_2xy
)+sumr
);
499 value
.i
= -(f
*sin_2xy
+sumi
);
500 value
= gem_cxmul(value
,gem_cmplx(exp_y2
*cos_2xy
,exp_y2
*sin_2xy
));
504 /* W(x,z) exp(-x^2)*omega(x+z) */
505 static gem_complex
gem_cW(double x
, gem_complex z
){ return(gem_cxrmul(gem_comega(gem_cxradd(z
,x
)),exp(-x
*x
))); }
507 /* ------------ end of [cr]error.c ------------ */
509 /*_ REVERSIBLE GEMINATE RECOMBINATION
511 * Here are the functions for reversible geminate recombination. */
513 /* Changes the unit from square cm per s to square Ångström per ps,
514 * since Omers code uses the latter units while g_mds outputs the former.
515 * g_hbond expects a diffusion coefficent given in square cm per s. */
516 static double sqcm_per_s_to_sqA_per_ps (real D
) {
517 fprintf(stdout
, "Diffusion coefficient is %f A^2/ps\n", D
*1e4
);
518 return (double)(D
*1e4
);
522 static double eq10v2(double theoryCt
[], double time
[], int manytimes
,
523 double ka
, double kd
, t_gemParams
*params
)
525 /* Finding the 3 roots */
527 double kD
, D
, r
, a
, b
, c
, tsqrt
, sumimaginary
;
532 part1
, part2
, part3
, part4
;
537 a
= (1.0 + ka
/kD
) * sqrt(D
)/r
;
541 gem_solve(&alpha
, &beta
, &gamma
, a
, b
, c
);
542 /* Finding the 3 roots */
545 part1
= gem_cxmul(alpha
, gem_cxmul(gem_cxadd(beta
, gamma
), gem_cxsub(beta
, gamma
))); /* 1(2+3)(2-3) */
546 part2
= gem_cxmul(beta
, gem_cxmul(gem_cxadd(gamma
, alpha
), gem_cxsub(gamma
, alpha
))); /* 2(3+1)(3-1) */
547 part3
= gem_cxmul(gamma
, gem_cxmul(gem_cxadd(alpha
, beta
) , gem_cxsub(alpha
, beta
))); /* 3(1+2)(1-2) */
548 part4
= gem_cxmul(gem_cxsub(gamma
, alpha
), gem_cxmul(gem_cxsub(alpha
, beta
), gem_cxsub(beta
, gamma
))); /* (3-1)(1-2)(2-3) */
551 #pragma omp parallel for \
552 private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4), \
553 reduction(+:sumimaginary), \
557 for (i
=0; i
<manytimes
; i
++){
558 tsqrt
= sqrt(time
[i
]);
559 oma
= gem_comega(gem_cxrmul(alpha
, tsqrt
));
560 c1
= gem_cxmul(oma
, gem_cxdiv(part1
, part4
));
561 omb
= gem_comega(gem_cxrmul(beta
, tsqrt
));
562 c2
= gem_cxmul(omb
, gem_cxdiv(part2
, part4
));
563 omc
= gem_comega(gem_cxrmul(gamma
, tsqrt
));
564 c3
= gem_cxmul(omc
, gem_cxdiv(part3
, part4
));
565 c4
.r
= c1
.r
+c2
.r
+c3
.r
;
566 c4
.i
= c1
.i
+c2
.i
+c3
.i
;
568 sumimaginary
+= c4
.i
* c4
.i
;
575 /* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
576 static double getLogIndex(const int i
, const t_gemParams
*params
)
578 return (exp(((double)(i
)) * params
->logQuota
) -1);
581 extern t_gemParams
*init_gemParams(const double sigma
, const double D
,
582 const real
*t
, const int len
, const int nFitPoints
,
583 const real begFit
, const real endFit
,
584 const real ballistic
, const int nBalExp
, const bool bDt
)
591 /* A few hardcoded things here. For now anyway. */
593 /* p->ka_max = 100; */
601 /* p->lifetime = 0; */
602 p
->sigma
= sigma
*10; /* Omer uses Ã…, not nm */
603 /* p->lsq_old = 99999; */
604 p
->D
= sqcm_per_s_to_sqA_per_ps(D
);
605 p
->kD
= 4*acos(-1.0)*sigma
*p
->D
;
608 /* Parameters used by calcsquare(). Better to calculate them
609 * here than in calcsquare every time it's called. */
611 /* p->logAfterTime = logAfterTime; */
612 tDelta
= (t
[len
-1]-t
[0]) / len
;
614 gmx_fatal(FARGS
, "Time between frames is non-positive!");
619 p
->nExpFit
= nBalExp
;
620 /* p->nLin = logAfterTime / tDelta; */
621 p
->nFitPoints
= nFitPoints
;
624 p
->logQuota
= (double)(log(p
->len
))/(p
->nFitPoints
-1);
625 /* if (p->nLin <= 0) { */
626 /* fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
630 /* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
631 /* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
632 /* p->logPF = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
633 /* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */
635 /* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
636 p
->ballistic
= ballistic
;
641 /* There was a misunderstanding regarding the fitting. From our
642 * recent correspondence it appears that Omer's code require
643 * the ACF data on a log-scale and does not operate on the raw data.
644 * This needs to be redone in gemFunc_residual() as well as in the
645 * t_gemParams structure. */
647 static double gemFunc_residual2(const gsl_vector
*p
, void *data
)
650 int i
, iLog
, nLin
, nFitPoints
, nData
;
651 double r
, residual2
, *ctTheory
, *y
;
653 GD
= (gemFitData
*)data
;
654 nLin
= GD
->params
->nLin
;
655 nFitPoints
= GD
->params
->nFitPoints
;
658 ctTheory
= GD
->ctTheory
;
661 /* Now, we'd save a lot of time by not calculating eq10v2 for all
662 * time points, but only those that we sample to calculate the mse.
665 eq10v2(GD
->ctTheory
, GD
->doubleLogTime
/* GD->time */, nFitPoints
/* GD->nData */,
666 gsl_vector_get(p
, 0), gsl_vector_get(p
, 1),
669 /* Removing a bunch of points from the log-part. */
671 #pragma omp parallel for schedule(dynamic) \
672 firstprivate(nData, ctTheory, y, nFitPoints) \
673 private (i, iLog, r) \
674 reduction(+:residual2) \
677 for(i
=0; i
<nFitPoints
; i
++)
679 iLog
= GD
->logtime
[i
];
680 r
= log(ctTheory
[i
/* iLog */]);
681 residual2
+= sqr(r
-log(y
[iLog
]));
683 residual2
/= nFitPoints
; /* Not really necessary for the fitting, but gives more meaning to the returned data. */
684 /* printf("ka=%3.5f\tkd=%3.5f\trmse=%3.5f\n", gsl_vector_get(p,0), gsl_vector_get(p,1), residual2); */
687 /* for (i=0; i<nLin; i++) { */
688 /* /\* Linear part ----------*\/ */
689 /* r = ctTheory[i]; */
690 /* residual2 += sqr(r-y[i]); */
691 /* /\* Log part -------------*\/ */
692 /* iLog = GETLOGINDEX(i, GD->params); */
693 /* if (iLog >= nData) */
694 /* gmx_fatal(FARGS, "log index out of bounds: %i", iLog); */
695 /* r = ctTheory[iLog]; */
696 /* residual2 += sqr(r-y[iLog]); */
699 /* residual2 /= GD->n; */
700 /* return residual2; */
704 static real
* d2r(const double *d
, const int nn
);
706 extern real
fitGemRecomb(double *ct
, double *time
, double **ctFit
,
707 const int nData
, t_gemParams
*params
)
710 int nThreads
, i
, iter
, status
, maxiter
;
711 real size
, d2
, tol
, *dumpdata
;
714 char *dumpstr
, dumpname
[128];
716 /* nmsimplex2 had convergence problems prior to gsl v1.14,
717 * but it's O(N) instead of O(N) operations, so let's use it if v >= 1.14 */
719 gsl_multimin_fminimizer
*s
;
720 gsl_vector
*x
,*dx
; /* parameters and initial step size */
721 gsl_multimin_function fitFunc
;
722 #ifdef GSL_MAJOR_VERSION
723 #ifdef GSL_MINOR_VERSION
724 #if ((GSL_MAJOR_VERSION == 1 && GSL_MINOR_VERSION >= 14) || \
725 (GSL_MAJOR_VERSION > 1))
726 const gsl_multimin_fminimizer_type
*T
= gsl_multimin_fminimizer_nmsimplex2
;
728 const gsl_multimin_fminimizer_type
*T
= gsl_multimin_fminimizer_nmsimplex
;
730 #endif /* GSL_MINOR_VERSION */
732 const gsl_multimin_fminimizer_type
*T
= gsl_multimin_fminimizer_nmsimplex
;
733 #endif /* GSL_MAJOR_VERSION */
734 fprintf(stdout
, "Will fit ka and kd to the ACF according to the reversible geminate recombination model.\n");
735 #else /* HAVE_LIBGSL */
736 fprintf(stderr
, "Sorry, can't do reversible geminate recombination without gsl. "
737 "Recompile using --with-gsl.\n");
739 #endif /* HAVE_LIBGSL */
743 nThreads
= omp_get_num_procs();
744 omp_set_num_threads(nThreads
);
745 fprintf(stdout
, "We will be using %i threads.\n", nThreads
);
753 p
= 2; /* Number of parameters to fit. ka and kd. */
754 n
= params
->nFitPoints
; /* params->nLin*2 */; /* Number of points in the reduced dataset */
758 fprintf(stderr
, "Fitting of D is not implemented yet. It must be provided on the command line.\n");
763 /* fprintf(stderr, "Reduced data set larger than the complete data set!\n"); */
766 snew(dumpdata
, nData
);
772 snew(GD
->ctTheory
, nData
);
778 GD
->tDelta
= time
[1]-time
[0];
781 snew(GD
->logtime
,params
->nFitPoints
);
782 snew(GD
->doubleLogTime
,params
->nFitPoints
);
784 for (i
=0; i
<params
->nFitPoints
; i
++)
786 GD
->doubleLogTime
[i
] = (double)(getLogIndex(i
, params
));
787 GD
->logtime
[i
] = (int)(GD
->doubleLogTime
[i
]);
788 GD
->doubleLogTime
[i
]*=GD
->tDelta
;
790 if (GD
->logtime
[i
] >= nData
)
792 fprintf(stderr
, "Ayay. It seems we're indexing out of bounds.\n");
793 params
->nFitPoints
= i
;
797 fitFunc
.f
= &gemFunc_residual2
;
799 fitFunc
.params
= (void*)GD
;
801 x
= gsl_vector_alloc (fitFunc
.n
);
802 dx
= gsl_vector_alloc (fitFunc
.n
);
803 gsl_vector_set (x
, 0, 25);
804 gsl_vector_set (x
, 1, 0.5);
805 gsl_vector_set (dx
, 0, 0.1);
806 gsl_vector_set (dx
, 1, 0.01);
809 s
= gsl_multimin_fminimizer_alloc (T
, fitFunc
.n
);
810 gsl_multimin_fminimizer_set (s
, &fitFunc
, x
, dx
);
812 gsl_vector_free (dx
);
816 status
= gsl_multimin_fminimizer_iterate (s
);
819 gmx_fatal(FARGS
,"Something went wrong in the iteration in minimizer %s:\n \"%s\"\n",
820 gsl_multimin_fminimizer_name(s
), gsl_strerror(status
));
822 d2
= gsl_multimin_fminimizer_minimum(s
);
823 size
= gsl_multimin_fminimizer_size(s
);
824 params
->ka
= gsl_vector_get (s
->x
, 0);
825 params
->kd
= gsl_vector_get (s
->x
, 1);
829 fprintf(stderr
, "%s\n", gsl_strerror(status
));
833 status
= gsl_multimin_test_size(size
,tol
);
835 if (status
== GSL_SUCCESS
) {
836 fprintf(stdout
, "Converged to minimum at\n");
839 printf ("iter %5d: ka = %2.5f kd = %2.5f f() = %7.3f size = %.3f chi2 = %2.5f\n",
847 eq10v2(GD
->ctTheory
, time
, nData
, params
->ka
, params
->kd
, params
);
848 sprintf(dumpname
, "Iter_%i.xvg", iter
);
849 for(i
=0; i
<GD
->nData
; i
++)
850 dumpdata
[i
] = (real
)(GD
->ctTheory
[i
]);
851 dumpN(dumpdata
, GD
->nData
, dumpname
);
854 while ((status
== GSL_CONTINUE
) && (iter
< maxiter
));
856 /* /\* Calculate the theoretical ACF from the parameters one last time. *\/ */
857 eq10v2(GD
->ctTheory
, time
, nData
, params
->ka
, params
->kd
, params
);
858 *ctFit
= GD
->ctTheory
;
861 gsl_multimin_fminimizer_free (s
);
866 #endif /* HAVE_LIBGSL */
870 static int balFunc_f(const gsl_vector
*x
, void *data
, gsl_vector
*f
)
872 /* C + sum{ A_i * exp(-B_i * t) }*/
876 double *y
, *A
, *B
, C
, /* There are the parameters to be optimized. */
879 BD
= (balData
*)data
;
886 for (i
= 0; i
<nexp
; i
++)
888 A
[i
] = gsl_vector_get(x
, i
*2);
889 B
[i
] = gsl_vector_get(x
, i
*2+1);
891 C
= gsl_vector_get(x
, nexp
*2);
897 for (j
=0; j
<nexp
; j
++) {
898 ct
+= A
[j
] * exp(B
[j
] * t
);
901 gsl_vector_set (f
, i
, ct
- y
[i
]);
906 /* The derivative stored in jacobian form (J)*/
907 static int balFunc_df(const gsl_vector
*params
, void *data
, gsl_matrix
*J
)
911 double *y
, *A
, *B
, C
, /* There are the parameters. */
923 for (i
=0; i
<nexp
; i
++)
925 A
[i
] = gsl_vector_get(params
, i
*2);
926 B
[i
] = gsl_vector_get(params
, i
*2+1);
928 C
= gsl_vector_get(params
, nexp
*2);
932 for (j
=0; j
<nexp
; j
++)
934 gsl_matrix_set (J
, i
, j
*2, exp(B
[j
]*t
)); /* df(t)/dA_j */
935 gsl_matrix_set (J
, i
, j
*2+1, A
[j
]*t
*exp(B
[j
]*t
)); /* df(t)/dB_j */
937 gsl_matrix_set (J
, i
, nexp
*2, 1); /* df(t)/dC */
942 /* Calculation of the function and its derivative */
943 static int balFunc_fdf(const gsl_vector
*params
, void *data
,
944 gsl_vector
*f
, gsl_matrix
*J
)
946 balFunc_f(params
, data
, f
);
947 balFunc_df(params
, data
, J
);
950 #endif /* HAVE_LIBGSL */
952 /* Removes the ballistic term from the beginning of the ACF,
953 * just like in Omer's paper.
955 extern void takeAwayBallistic(double *ct
, double *t
, int len
, real tMax
, int nexp
, bool bDerivative
)
958 /* Use nonlinear regression with GSL instead.
959 * Fit with 4 exponentials and one constant term,
960 * subtract the fatest exponential. */
962 int nData
,i
,status
, iter
;
964 double *guess
, /* Initial guess. */
965 *A
, /* The fitted parameters. (A1, B1, A2, B2,... C) */
975 } while (t
[nData
]<tMax
+t
[0] && nData
<len
);
977 p
= nexp
*2+1; /* Number of parameters. */
980 const gsl_multifit_fdfsolver_type
*T
981 = gsl_multifit_fdfsolver_lmsder
;
983 gsl_multifit_fdfsolver
*s
; /* The solver itself. */
984 gsl_multifit_function_fdf fitFunction
; /* The function to be fitted. */
985 gsl_matrix
*covar
; /* Covariance matrix for the parameters.
986 * We'll not use the result, though. */
987 gsl_vector_view theParams
;
992 } while (t
[nData
]<tMax
+t
[0] && nData
<len
);
999 covar
= gsl_matrix_alloc (p
, p
);
1001 /* Set up an initial gess for the parameters.
1002 * The solver is somewhat sensitive to the initial guess,
1003 * but this worked fine for a TIP3P box with -geminate dd
1004 * EDIT: In fact, this seems like a good starting pont for other watermodels too. */
1005 for (i
=0; i
<nexp
; i
++)
1008 guess
[i
*2+1] = -0.5 + (((double)i
)/nexp
- 0.5)*0.3;
1010 guess
[nexp
* 2] = 0.01;
1012 theParams
= gsl_vector_view_array(guess
, p
);
1017 BD
->tDelta
= t
[1]-t
[0];
1020 fitFunction
.f
= &balFunc_f
;
1021 fitFunction
.df
= &balFunc_df
;
1022 fitFunction
.fdf
= &balFunc_fdf
;
1023 fitFunction
.n
= nData
;
1025 fitFunction
.params
= BD
;
1027 s
= gsl_multifit_fdfsolver_alloc (T
, nData
, p
);
1029 gmx_fatal(FARGS
, "Could not set up the nonlinear solver.");
1031 gsl_multifit_fdfsolver_set(s
, &fitFunction
, &theParams
.vector
);
1033 /* \=============================================/ */
1039 status
= gsl_multifit_fdfsolver_iterate (s
);
1043 status
= gsl_multifit_test_delta (s
->dx
, s
->x
, 1e-4, 1e-4);
1045 while (iter
< 5000 && status
== GSL_CONTINUE
);
1049 fprintf(stderr
, "The non-linear fitting did not converge in 5000 steps.\n"
1050 "Check the quality of the fit!\n");
1054 fprintf(stderr
, "Non-linear fitting of ballistic term converged in %d steps.\n\n", (int)iter
);
1056 for (i
=0; i
<nexp
; i
++) {
1057 fprintf(stdout
, "%c * exp(%c * t) + ", 'A'+(char)i
*2, 'B'+(char)i
*2);
1060 fprintf(stdout
, "%c\n", 'A'+(char)nexp
*2);
1061 fprintf(stdout
, "Here are the actual numbers for A-%c:\n", 'A'+nexp
*2);
1063 for (i
=0; i
<nexp
; i
++)
1065 A
[i
*2] = gsl_vector_get(s
->x
, i
*2);
1066 A
[i
*2+1] = gsl_vector_get(s
->x
, i
*2+1);
1067 fprintf(stdout
, " %g*exp(%g * x) +", A
[i
*2], A
[i
*2+1]);
1069 A
[i
*2] = gsl_vector_get(s
->x
, i
*2); /* The last and constant term */
1070 fprintf(stdout
, " %g\n", A
[i
*2]);
1074 /* Implement some check for parameter quality */
1075 for (i
=0; i
<nexp
; i
++)
1077 if (A
[i
*2]<0 || A
[i
*2]>1) {
1078 fprintf(stderr
, "WARNING: ----------------------------------\n"
1079 " | A coefficient does not lie within [0,1].\n"
1080 " | This may or may not be a problem.\n"
1081 " | Double check the quality of the fit!\n");
1084 fprintf(stderr
, "WARNING: ----------------------------------\n"
1085 " | One factor in the exponent is positive.\n"
1086 " | This could be a problem if the coefficient\n"
1087 " | is large. Double check the quality of the fit!\n");
1090 if (A
[i
*2]<0 || A
[i
*2]>1) {
1091 fprintf(stderr
, "WARNING: ----------------------------------\n"
1092 " | The constant term does not lie within [0,1].\n"
1093 " | This may or may not be a problem.\n"
1094 " | Double check the quality of the fit!\n");
1097 /* Sort the terms */
1098 sorted
= (nexp
> 1) ? FALSE
: TRUE
;
1102 for (i
=0;i
<nexp
-1;i
++)
1104 ddt
[0] = A
[i
*2] * A
[i
*2+1];
1105 ddt
[1] =A
[i
*2+2] * A
[i
*2+3];
1107 if ((bDerivative
&& (ddt
[0]<0 && ddt
[1]<0 && ddt
[0]>ddt
[1])) || /* Compare derivative at t=0... */
1108 (!bDerivative
&& (A
[i
*2+1] > A
[i
*2+3]))) /* Or just the coefficient in the exponent */
1111 a
[0] = A
[i
*2]; /* coefficient */
1112 a
[1] = A
[i
*2+1]; /* parameter in the exponent */
1115 A
[i
*2+1] = A
[i
*2+3];
1123 /* Subtract the fastest component */
1124 fprintf(stdout
, "Fastest component is %g * exp(%g * t)\n"
1125 "Subtracting fastest component from ACF.\n", A
[0], A
[1]);
1127 for (i
=0; i
<len
; i
++) {
1128 ct
[i
] = (ct
[i
] - A
[0] * exp(A
[1] * i
*BD
->tDelta
)) / (1-A
[0]);
1134 gsl_multifit_fdfsolver_free(s
);
1135 gsl_matrix_free(covar
);
1139 /* We have no gsl. */
1140 fprintf(stderr
, "Sorry, can't take away ballistic component without gsl. "
1141 "Recompile using --with-gsl.\n");
1143 #endif /* HAVE_LIBGSL */
1147 extern void dumpN(const real
*e
, const int nn
, char *fn
)
1149 /* For debugging only */
1152 char standardName
[] = "Nt.xvg";
1160 "@ xaxis label \"Frame\"\n"
1161 "@ yaxis label \"N\"\n"
1162 "@ s0 line type 3\n");
1164 for (i
=0; i
<nn
; i
++) {
1165 fprintf(f
, "%-10i %-g\n", i
, e
[i
]);
1171 static real
* d2r(const double *d
, const int nn
)
1177 for (i
=0; i
<nn
; i
++)