src/gromacs/gmxana/geminate.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 4.5
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #ifdef HAVE_LIBGSL
  40 #include <gsl/gsl_rng.h>
  41 #include <gsl/gsl_randist.h>
  42 #include <gsl/gsl_vector.h>
  43 #include <gsl/gsl_blas.h>
  44 #include <gsl/gsl_multimin.h>
  45 #include <gsl/gsl_multifit_nlin.h>
  46 #include <gsl/gsl_sf.h>
  47 #include <gsl/gsl_version.h>
  48 #endif
  49
  50 #include <math.h>
  51 #include "typedefs.h"
  52 #include "smalloc.h"
  53 #include "vec.h"
  54 #include "geminate.h"
  55 #include "gmx_omp.h"
  56
  57
  58 /* The first few sections of this file contain functions that were adopted,
  59  * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
  60  * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
  61  * This is also the case with the function eq10v2().
  62  *
  63  * The parts menetioned in the previous paragraph were contributed under the BSD license.
  64  */
  65
  66
  67 /* This first part is from complex.c which I recieved from Omer Markowitch.
  68  * - Erik Marklund
  69  *
  70  * ------------- from complex.c ------------- */
  71
  72 /* Complexation of a paired number (r,i)                                     */
  73 static gem_complex gem_cmplx(double r, double i)
  74 {
  75     gem_complex value;
  76     value.r = r;
  77     value.i = i;
  78     return value;
  79 }
  80
  81 /* Complexation of a real number, x */
  82 static gem_complex gem_c(double x)
  83 {
  84     gem_complex value;
  85     value.r = x;
  86     value.i = 0;
  87     return value;
  88 }
  89
  90 /* Magnitude of a complex number z                                           */
  91 static double gem_cx_abs(gem_complex z)
  92 {
  93     return (sqrt(z.r*z.r+z.i*z.i));
  94 }
  95
  96 /* Addition of two complex numbers z1 and z2                                 */
  97 static gem_complex gem_cxadd(gem_complex z1, gem_complex z2)
  98 {
  99     gem_complex value;
 100     value.r = z1.r+z2.r;
 101     value.i = z1.i+z2.i;
 102     return value;
 103 }
 104
 105 /* Addition of a complex number z1 and a real number r */
 106 static gem_complex gem_cxradd(gem_complex z, double r)
 107 {
 108     gem_complex value;
 109     value.r = z.r + r;
 110     value.i = z.i;
 111     return value;
 112 }
 113
 114 /* Subtraction of two complex numbers z1 and z2                              */
 115 static gem_complex gem_cxsub(gem_complex z1, gem_complex z2)
 116 {
 117     gem_complex value;
 118     value.r = z1.r-z2.r;
 119     value.i = z1.i-z2.i;
 120     return value;
 121 }
 122
 123 /* Multiplication of two complex numbers z1 and z2                           */
 124 static gem_complex gem_cxmul(gem_complex z1, gem_complex z2)
 125 {
 126     gem_complex value;
 127     value.r = z1.r*z2.r-z1.i*z2.i;
 128     value.i = z1.r*z2.i+z1.i*z2.r;
 129     return value;
 130 }
 131
 132 /* Square of a complex number z                                              */
 133 static gem_complex gem_cxsq(gem_complex z)
 134 {
 135     gem_complex value;
 136     value.r = z.r*z.r-z.i*z.i;
 137     value.i = z.r*z.i*2.;
 138     return value;
 139 }
 140
 141 /* multiplication of a complex number z and a real number r */
 142 static gem_complex gem_cxrmul(gem_complex z, double r)
 143 {
 144     gem_complex value;
 145     value.r = z.r*r;
 146     value.i = z.i*r;
 147     return value;
 148 }
 149
 150 /* Division of two complex numbers z1 and z2                                 */
 151 static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2)
 152 {
 153     gem_complex value;
 154     double      num;
 155     num = z2.r*z2.r+z2.i*z2.i;
 156     if (num == 0.)
 157     {
 158         fprintf(stderr, "ERROR in gem_cxdiv function\n");
 159     }
 160     value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num;
 161     return value;
 162 }
 163
 164 /* division of a complex z number by a real number x */
 165 static gem_complex gem_cxrdiv(gem_complex z, double r)
 166 {
 167     gem_complex value;
 168     value.r = z.r/r;
 169     value.i = z.i/r;
 170     return value;
 171 }
 172
 173 /* division of a real number r by a complex number x */
 174 static gem_complex gem_rxcdiv(double r, gem_complex z)
 175 {
 176     gem_complex value;
 177     double      f;
 178     f       = r/(z.r*z.r+z.i*z.i);
 179     value.r = f*z.r;
 180     value.i = -f*z.i;
 181     return value;
 182 }
 183
 184 /* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
 185 static gem_complex gem_cxdexp(gem_complex z)
 186 {
 187     gem_complex value;
 188     double      exp_z_r;
 189     exp_z_r = exp(z.r);
 190     value.r = exp_z_r*cos(z.i);
 191     value.i = exp_z_r*sin(z.i);
 192     return value;
 193 }
 194
 195 /* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z),                  */
 196 /*                                  where -PI < Arg(z) < PI                  */
 197 static gem_complex gem_cxlog(gem_complex z)
 198 {
 199     gem_complex value;
 200     double      mag2;
 201     mag2 = z.r*z.r+z.i*z.i;
 202     if (mag2 < 0.)
 203     {
 204         fprintf(stderr, "ERROR in gem_cxlog func\n");
 205     }
 206     value.r = log(sqrt(mag2));
 207     if (z.r == 0.)
 208     {
 209         value.i = PI/2.;
 210         if (z.i < 0.)
 211         {
 212             value.i = -value.i;
 213         }
 214     }
 215     else
 216     {
 217         value.i = atan2(z.i, z.r);
 218     }
 219     return value;
 220 }
 221
 222 /* Square root of a complex number z = |z| exp(I*the) -- z^(1/2)             */
 223 /*                               z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */
 224 /*                               where 0 < the < 2*PI                        */
 225 static gem_complex gem_cxdsqrt(gem_complex z)
 226 {
 227     gem_complex value;
 228     double      sq;
 229     sq      = gem_cx_abs(z);
 230     value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */
 231     value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */
 232     if (z.i < 0.)
 233     {
 234         value.r = -value.r;
 235     }
 236     return value;
 237 }
 238
 239 /* Complex power of a complex number z1^z2                                   */
 240 static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2)
 241 {
 242     gem_complex value;
 243     value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2));
 244     return value;
 245 }
 246
 247 /* ------------ end of complex.c ------------ */
 248
 249 /* This next part was derived from cubic.c, also received from Omer Markovitch.
 250  * ------------- from cubic.c ------------- */
 251
 252 /* Solver for a cubic equation: x^3-a*x^2+b*x-c=0                            */
 253 static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam,
 254                       double a, double b, double c)
 255 {
 256     double      t1, t2, two_3, temp;
 257     gem_complex ctemp, ct3;
 258
 259     two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c;
 260     temp  = 4.*t1*t1*t1+t2*t2;
 261
 262     ctemp = gem_cmplx(temp, 0.);   ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp));
 263     ct3   = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.));
 264
 265     ctemp = gem_rxcdiv(-two_3*t1/3., ct3);
 266     ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3));
 267
 268     *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp);
 269
 270     ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.))));
 271     ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam));
 272     ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam);
 273     *al   = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
 274     *be   = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
 275 }
 276
 277 /* ------------ end of cubic.c ------------ */
 278
 279 /* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
 280  * ------------- from [cr]error.c ------------- */
 281
 282 /************************************************************/
 283 /* Real valued error function and related functions         */
 284 /* x, y     : real variables                                */
 285 /* erf(x)   : error function                                */
 286 /* erfc(x)  : complementary error function                  */
 287 /* omega(x) : exp(x*x)*erfc(x)                              */
 288 /* W(x,y)   : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
 289 /************************************************************/
 290
 291 /*---------------------------------------------------------------------------*/
 292 /* Utilzed the series approximation of erf(x)                                */
 293 /* Relative error=|err(x)|/erf(x)<EPS                                        */
 294 /* Handbook of Mathematical functions, Abramowitz, p 297                     */
 295 /* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
 296 /*---------------------------------------------------------------------------*/
 297 /*  This gives erfc(z) correctly only upto >10-15 */
 298
 299 static double gem_erf(double x)
 300 {
 301     double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12;
 302     temp = x;
 303     sum  = temp;
 304     xm   = 26.;
 305     x2   = x*x;
 306     x4   = x2*x2;
 307     x6   = x4*x2;
 308     x8   = x6*x2;
 309     x10  = x8*x2;
 310     x12  = x10*x2;
 311     exp2 = exp(-x2);
 312     if (x <= xm)
 313     {
 314         for (n = 1.; n <= 2000.; n += 1.)
 315         {
 316             temp *= 2.*x2/(2.*n+1.);
 317             sum  += temp;
 318             if (fabs(temp/sum) < 1.E-16)
 319             {
 320                 break;
 321             }
 322         }
 323
 324         if (n >= 2000.)
 325         {
 326             fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n);
 327         }
 328         sum *= 2./sPI*exp2;
 329     }
 330     else
 331     {
 332         /* from the asymptotic expansion of experfc(x) */
 333         sum = (1. - 0.5/x2 + 0.75/x4
 334                - 1.875/x6 + 6.5625/x8
 335                - 29.53125/x10 + 162.421875/x12)
 336             / sPI/x;
 337         sum *= exp2; /* now sum is erfc(x) */
 338         sum  = -sum+1.;
 339     }
 340     return x >= 0.0 ? sum : -sum;
 341 }
 342
 343 /* Result --> Alex's code for erfc and experfc behaves better than this      */
 344 /* complementray error function.                Returns 1.-erf(x)            */
 345 static double gem_erfc(double x)
 346 {
 347     double t, z, ans;
 348     z = fabs(x);
 349     t = 1.0/(1.0+0.5*z);
 350
 351     ans = t * exp(-z*z - 1.26551223 +
 352                   t*(1.00002368 +
 353                      t*(0.37409196 +
 354                         t*(0.09678418 +
 355                            t*(-0.18628806 +
 356                               t*(0.27886807 +
 357                                  t*(-1.13520398 +
 358                                     t*(1.48851587 +
 359                                        t*(-0.82215223 +
 360                                           t*0.17087277)))))))));
 361
 362     return x >= 0.0 ? ans : 2.0-ans;
 363 }
 364
 365 /* omega(x)=exp(x^2)erfc(x)                                                  */
 366 static double gem_omega(double x)
 367 {
 368     double xm, ans, xx, x4, x6, x8, x10, x12;
 369     xm  = 26;
 370     xx  = x*x;
 371     x4  = xx*xx;
 372     x6  = x4*xx;
 373     x8  = x6*xx;
 374     x10 = x8*xx;
 375     x12 = x10*xx;
 376
 377     if (x <= xm)
 378     {
 379         ans = exp(xx)*gem_erfc(x);
 380     }
 381     else
 382     {
 383         /* Asymptotic expansion */
 384         ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x;
 385     }
 386     return ans;
 387 }
 388
 389 /*---------------------------------------------------------------------------*/
 390 /* Utilzed the series approximation of erf(z=x+iy)                           */
 391 /* Relative error=|err(z)|/|erf(z)|<EPS                                      */
 392 /* Handbook of Mathematical functions, Abramowitz, p 299                     */
 393 /* comega(z=x+iy)=cexp(z^2)*cerfc(z)                                         */
 394 /*---------------------------------------------------------------------------*/
 395 static gem_complex gem_comega(gem_complex z)
 396 {
 397     gem_complex value;
 398     double      x, y;
 399     double      sumr, sumi, n, n2, f, temp, temp1;
 400     double      x2, y2, cos_2xy, sin_2xy, cosh_2xy, sinh_2xy, cosh_ny, sinh_ny, exp_y2;
 401
 402     x        = z.r;
 403     y        = z.i;
 404     x2       = x*x;
 405     y2       = y*y;
 406     sumr     = 0.;
 407     sumi     = 0.;
 408     cos_2xy  = cos(2.*x*y);
 409     sin_2xy  = sin(2.*x*y);
 410     cosh_2xy = cosh(2.*x*y);
 411     sinh_2xy = sinh(2.*x*y);
 412     exp_y2   = exp(-y2);
 413
 414     for (n = 1.0, temp = 0.; n <= 2000.; n += 1.0)
 415     {
 416         n2      = n*n;
 417         cosh_ny = cosh(n*y);
 418         sinh_ny = sinh(n*y);
 419         f       = exp(-n2/4.)/(n2+4.*x2);
 420         /* if(f<1.E-200) break; */
 421         sumr += (2.*x - 2.*x*cosh_ny*cos_2xy + n*sinh_ny*sin_2xy)*f;
 422         sumi += (2.*x*cosh_ny*sin_2xy + n*sinh_ny*cos_2xy)*f;
 423         temp1 = sqrt(sumr*sumr+sumi*sumi);
 424         if (fabs((temp1-temp)/temp1) < 1.E-16)
 425         {
 426             break;
 427         }
 428         temp = temp1;
 429     }
 430     if (n == 2000.)
 431     {
 432         fprintf(stderr, "iteration exceeds %lg\n", n);
 433     }
 434     sumr *= 2./PI;
 435     sumi *= 2./PI;
 436
 437     if (x != 0.)
 438     {
 439         f = 1./2./PI/x;
 440     }
 441     else
 442     {
 443         f = 0.;
 444     }
 445     value.r = gem_omega(x)-(f*(1.-cos_2xy)+sumr);
 446     value.i = -(f*sin_2xy+sumi);
 447     value   = gem_cxmul(value, gem_cmplx(exp_y2*cos_2xy, exp_y2*sin_2xy));
 448     return (value);
 449 }
 450
 451 /* ------------ end of [cr]error.c ------------ */
 452
 453 /*_ REVERSIBLE GEMINATE RECOMBINATION
 454  *
 455  * Here are the functions for reversible geminate recombination. */
 456
 457 /* Changes the unit from square cm per s to square Ångström per ps,
 458  * since Omers code uses the latter units while g_mds outputs the former.
 459  * g_hbond expects a diffusion coefficent given in square cm per s. */
 460 static double sqcm_per_s_to_sqA_per_ps (real D)
 461 {
 462     fprintf(stdout, "Diffusion coefficient is %f A^2/ps\n", D*1e4);
 463     return (double)(D*1e4);
 464 }
 465
 466
 467 static double eq10v2(double theoryCt[], double time[], int manytimes,
 468                      double ka, double kd, t_gemParams *params)
 469 {
 470     /* Finding the 3 roots */
 471     int    i;
 472     double kD, D, r, a, b, c, tsqrt, sumimaginary;
 473     gem_complex
 474            alpha, beta, gamma,
 475         c1, c2, c3, c4,
 476         oma, omb, omc,
 477         part1, part2, part3, part4;
 478
 479     kD = params->kD;
 480     D  = params->D;
 481     r  = params->sigma;
 482     a  = (1.0 + ka/kD) * sqrt(D)/r;
 483     b  = kd;
 484     c  = kd * sqrt(D)/r;
 485
 486     gem_solve(&alpha, &beta, &gamma, a, b, c);
 487     /* Finding the 3 roots */
 488
 489     sumimaginary = 0;
 490     part1        = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta,  gamma), gem_cxsub(beta,  gamma)));                 /* 1(2+3)(2-3) */
 491     part2        = gem_cxmul(beta,  gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha)));                 /* 2(3+1)(3-1) */
 492     part3        = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta)));                   /* 3(1+2)(1-2) */
 493     part4        = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */
 494
 495 #pragma omp parallel for \
 496     private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
 497     reduction(+:sumimaginary) \
 498     default(shared) \
 499     schedule(guided)
 500     for (i = 0; i < manytimes; i++)
 501     {
 502         tsqrt         = sqrt(time[i]);
 503         oma           = gem_comega(gem_cxrmul(alpha, tsqrt));
 504         c1            = gem_cxmul(oma, gem_cxdiv(part1, part4));
 505         omb           = gem_comega(gem_cxrmul(beta, tsqrt));
 506         c2            = gem_cxmul(omb, gem_cxdiv(part2, part4));
 507         omc           = gem_comega(gem_cxrmul(gamma, tsqrt));
 508         c3            = gem_cxmul(omc, gem_cxdiv(part3, part4));
 509         c4.r          = c1.r+c2.r+c3.r;
 510         c4.i          = c1.i+c2.i+c3.i;
 511         theoryCt[i]   = c4.r;
 512         sumimaginary += c4.i * c4.i;
 513     }
 514
 515     return sumimaginary;
 516
 517 } /* eq10v2 */
 518
 519 /* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
 520 static double getLogIndex(const int i, const t_gemParams *params)
 521 {
 522     return (exp(((double)(i)) * params->logQuota) -1);
 523 }
 524
 525 extern t_gemParams *init_gemParams(const double sigma, const double D,
 526                                    const real *t, const int len, const int nFitPoints,
 527                                    const real begFit, const real endFit,
 528                                    const real ballistic, const int nBalExp, const gmx_bool bDt)
 529 {
 530     double       tDelta;
 531     t_gemParams *p;
 532
 533     snew(p, 1);
 534
 535     /* A few hardcoded things here. For now anyway. */
 536 /*   p->ka_min   = 0; */
 537 /*   p->ka_max   = 100; */
 538 /*   p->dka      = 10; */
 539 /*   p->kd_min   = 0; */
 540 /*   p->kd_max   = 2; */
 541 /*   p->dkd      = 0.2; */
 542     p->ka       = 0;
 543     p->kd       = 0;
 544 /*   p->lsq      = -1; */
 545 /*   p->lifetime = 0; */
 546     p->sigma    = sigma*10; /* Omer uses Å, not nm */
 547 /*   p->lsq_old  = 99999; */
 548     p->D        = sqcm_per_s_to_sqA_per_ps(D);
 549     p->kD       = 4*acos(-1.0)*sigma*p->D;
 550
 551
 552     /* Parameters used by calcsquare(). Better to calculate them
 553      * here than in calcsquare every time it's called. */
 554     p->len = len;
 555 /*   p->logAfterTime = logAfterTime; */
 556     tDelta       = (t[len-1]-t[0]) / len;
 557     if (tDelta <= 0)
 558     {
 559         gmx_fatal(FARGS, "Time between frames is non-positive!");
 560     }
 561     else
 562     {
 563         p->tDelta = tDelta;
 564     }
 565
 566     p->nExpFit      = nBalExp;
 567 /*   p->nLin         = logAfterTime / tDelta; */
 568     p->nFitPoints   = nFitPoints;
 569     p->begFit       = begFit;
 570     p->endFit       = endFit;
 571     p->logQuota     = (double)(log(p->len))/(p->nFitPoints-1);
 572 /*   if (p->nLin <= 0) { */
 573 /*     fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
 574 /*     sfree(p); */
 575 /*     return NULL; */
 576 /*   } */
 577 /* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
 578 /* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
 579 /*   p->logPF    = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
 580 /* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */
 581
 582 /*   p->logMult      = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
 583     p->ballistic    =  ballistic;
 584     return p;
 585 }
 586
 587 /* There was a misunderstanding regarding the fitting. From our
 588  * recent correspondence it appears that Omer's code require
 589  * the ACF data on a log-scale and does not operate on the raw data.
 590  * This needs to be redone in gemFunc_residual() as well as in the
 591  * t_gemParams structure. */
 592 #ifdef HAVE_LIBGSL
 593 static double gemFunc_residual2(const gsl_vector *p, void *data)
 594 {
 595     gemFitData *GD;
 596     int         i, iLog, nLin, nFitPoints, nData;
 597     double      r, residual2, *ctTheory, *y;
 598
 599     GD         = (gemFitData *)data;
 600     nLin       = GD->params->nLin;
 601     nFitPoints = GD->params->nFitPoints;
 602     nData      = GD->nData;
 603     residual2  = 0;
 604     ctTheory   = GD->ctTheory;
 605     y          = GD->y;
 606
 607     /* Now, we'd save a lot of time by not calculating eq10v2 for all
 608      * time points, but only those that we sample to calculate the mse.
 609      */
 610
 611     eq10v2(GD->ctTheory, GD->doubleLogTime /* GD->time */, nFitPoints /* GD->nData */,
 612            gsl_vector_get(p, 0), gsl_vector_get(p, 1),
 613            GD->params);
 614
 615     fixGemACF(GD->ctTheory, nFitPoints);
 616
 617     /* Removing a bunch of points from the log-part. */
 618 #pragma omp parallel for schedule(dynamic) \
 619     firstprivate(nData, ctTheory, y, nFitPoints) \
 620     private (i, iLog, r) \
 621     reduction(+:residual2) \
 622     default(shared)
 623     for (i = 0; i < nFitPoints; i++)
 624     {
 625         iLog       = GD->logtime[i];
 626         r          = log(ctTheory[i /* iLog */]);
 627         residual2 += sqr(r-log(y[iLog]));
 628     }
 629     residual2 /= nFitPoints; /* Not really necessary for the fitting, but gives more meaning to the returned data. */
 630     /* printf("ka=%3.5f\tkd=%3.5f\trmse=%3.5f\n", gsl_vector_get(p,0), gsl_vector_get(p,1), residual2); */
 631     return residual2;
 632
 633 /*   for (i=0; i<nLin; i++) { */
 634 /*     /\* Linear part ----------*\/ */
 635 /*     r = ctTheory[i]; */
 636 /*     residual2 += sqr(r-y[i]); */
 637 /*     /\* Log part -------------*\/ */
 638 /*     iLog = GETLOGINDEX(i, GD->params); */
 639 /*     if (iLog >= nData) */
 640 /*       gmx_fatal(FARGS, "log index out of bounds: %i", iLog); */
 641 /*     r = ctTheory[iLog]; */
 642 /*     residual2 += sqr(r-y[iLog]); */
 643
 644 /*   } */
 645 /*   residual2 /= GD->n; */
 646 /*   return residual2; */
 647 }
 648 #endif
 649
 650 static real* d2r(const double *d, const int nn);
 651
 652 extern real fitGemRecomb(double *ct, double *time, double **ctFit,
 653                          const int nData, t_gemParams *params)
 654 {
 655
 656     int         nThreads, i, iter, status, maxiter;
 657     real        size, d2, tol, *dumpdata;
 658     size_t      p, n;
 659     gemFitData *GD;
 660     char       *dumpstr, dumpname[128];
 661
 662     /* nmsimplex2 had convergence problems prior to gsl v1.14,
 663      * but it's O(N) instead of O(N) operations, so let's use it if v >= 1.14 */
 664 #ifdef HAVE_LIBGSL
 665     gsl_multimin_fminimizer *s;
 666     gsl_vector              *x, *dx; /* parameters and initial step size */
 667     gsl_multimin_function    fitFunc;
 668 #ifdef GSL_MAJOR_VERSION
 669 #ifdef GSL_MINOR_VERSION
 670 #if ((GSL_MAJOR_VERSION == 1 && GSL_MINOR_VERSION >= 14) || \
 671     (GSL_MAJOR_VERSION > 1))
 672     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
 673 #else
 674     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
 675 #endif /* #if ... */
 676 #endif /* GSL_MINOR_VERSION */
 677 #else
 678     const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
 679 #endif /* GSL_MAJOR_VERSION */
 680     fprintf(stdout, "Will fit ka and kd to the ACF according to the reversible geminate recombination model.\n");
 681 #else  /* HAVE_LIBGSL */
 682     fprintf(stderr, "Sorry, can't do reversible geminate recombination without gsl. "
 683             "Recompile using --with-gsl.\n");
 684     return -1;
 685 #endif /* HAVE_LIBGSL */
 686
 687 #ifdef HAVE_LIBGSL
 688 #ifdef GMX_OPENMP
 689     nThreads = gmx_omp_get_max_threads();
 690     gmx_omp_set_num_threads(nThreads);
 691     fprintf(stdout, "We will be using %i threads.\n", nThreads);
 692 #endif
 693
 694     iter    = 0;
 695     status  = 0;
 696     maxiter = 100;
 697     tol     = 1e-10;
 698
 699     p = 2;                                        /* Number of parameters to fit. ka and kd.  */
 700     n = params->nFitPoints; /* params->nLin*2 */; /* Number of points in the reduced dataset  */
 701
 702     if (params->D <= 0)
 703     {
 704         fprintf(stderr, "Fitting of D is not implemented yet. It must be provided on the command line.\n");
 705         return -1;
 706     }
 707
 708 /*   if (nData<n) { */
 709 /*     fprintf(stderr, "Reduced data set larger than the complete data set!\n"); */
 710 /*     n=nData; */
 711 /*   } */
 712     snew(dumpdata, nData);
 713     snew(GD, 1);
 714
 715     GD->n        = n;
 716     GD->y        = ct;
 717     GD->ctTheory = NULL;
 718     snew(GD->ctTheory, nData);
 719     GD->LinLog = NULL;
 720     snew(GD->LinLog, n);
 721     GD->time   = time;
 722     GD->ka     = 0;
 723     GD->kd     = 0;
 724     GD->tDelta = time[1]-time[0];
 725     GD->nData  = nData;
 726     GD->params = params;
 727     snew(GD->logtime, params->nFitPoints);
 728     snew(GD->doubleLogTime, params->nFitPoints);
 729
 730     for (i = 0; i < params->nFitPoints; i++)
 731     {
 732         GD->doubleLogTime[i]  = (double)(getLogIndex(i, params));
 733         GD->logtime[i]        = (int)(GD->doubleLogTime[i]);
 734         GD->doubleLogTime[i] *= GD->tDelta;
 735
 736         if (GD->logtime[i] >= nData)
 737         {
 738             fprintf(stderr, "Ayay. It seems we're indexing out of bounds.\n");
 739             params->nFitPoints = i;
 740         }
 741     }
 742
 743     fitFunc.f      = &gemFunc_residual2;
 744     fitFunc.n      = 2;
 745     fitFunc.params = (void*)GD;
 746
 747     x  = gsl_vector_alloc (fitFunc.n);
 748     dx = gsl_vector_alloc (fitFunc.n);
 749     gsl_vector_set (x,  0, 25);
 750     gsl_vector_set (x,  1, 0.5);
 751     gsl_vector_set (dx, 0, 0.1);
 752     gsl_vector_set (dx, 1, 0.01);
 753
 754
 755     s = gsl_multimin_fminimizer_alloc (T, fitFunc.n);
 756     gsl_multimin_fminimizer_set (s, &fitFunc, x, dx);
 757     gsl_vector_free (x);
 758     gsl_vector_free (dx);
 759
 760     do
 761     {
 762         iter++;
 763         status = gsl_multimin_fminimizer_iterate (s);
 764
 765         if (status != 0)
 766         {
 767             gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s:\n \"%s\"\n",
 768                       gsl_multimin_fminimizer_name(s), gsl_strerror(status));
 769         }
 770
 771         d2         = gsl_multimin_fminimizer_minimum(s);
 772         size       = gsl_multimin_fminimizer_size(s);
 773         params->ka = gsl_vector_get (s->x, 0);
 774         params->kd = gsl_vector_get (s->x, 1);
 775
 776         if (status)
 777         {
 778             fprintf(stderr, "%s\n", gsl_strerror(status));
 779             break;
 780         }
 781
 782         status = gsl_multimin_test_size(size, tol);
 783
 784         if (status == GSL_SUCCESS)
 785         {
 786             fprintf(stdout, "Converged to minimum at\n");
 787         }
 788
 789         printf ("iter %5d: ka = %2.5f  kd = %2.5f  f() = %7.3f  size = %.3f  chi2 = %2.5f\n",
 790                 iter,
 791                 params->ka,
 792                 params->kd,
 793                 s->fval, size, d2);
 794
 795         if (iter%1 == 0)
 796         {
 797             eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params);
 798             /* fixGemACF(GD->ctTheory, nFitPoints); */
 799             sprintf(dumpname, "Iter_%i.xvg", iter);
 800             for (i = 0; i < GD->nData; i++)
 801             {
 802                 dumpdata[i] = (real)(GD->ctTheory[i]);
 803                 if (!gmx_isfinite(dumpdata[i]))
 804                 {
 805                     gmx_fatal(FARGS, "Non-finite value in acf.");
 806                 }
 807             }
 808             dumpN(dumpdata, GD->nData, dumpname);
 809         }
 810     }
 811     while ((status == GSL_CONTINUE) && (iter < maxiter));
 812
 813     /*   /\* Calculate the theoretical ACF from the parameters one last time. *\/ */
 814     eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params);
 815     *ctFit = GD->ctTheory;
 816
 817     sfree(GD);
 818     gsl_multimin_fminimizer_free (s);
 819
 820
 821     return d2;
 822
 823 #endif /* HAVE_LIBGSL */
 824 }
 825
 826 #ifdef HAVE_LIBGSL
 827 static int balFunc_f(const gsl_vector *x, void *data, gsl_vector *f)
 828 {
 829     /* C + sum{ A_i * exp(-B_i * t) }*/
 830
 831     balData *BD;
 832     int      n, i, j, nexp;
 833     double  *y, *A, *B, C,     /* There are the parameters to be optimized. */
 834              t, ct;
 835
 836     BD   = (balData *)data;
 837     n    = BD->n;
 838     nexp = BD->nexp;
 839     y    = BD->y,
 840     snew(A, nexp);
 841     snew(B, nexp);
 842
 843     for (i = 0; i < nexp; i++)
 844     {
 845         A[i] = gsl_vector_get(x, i*2);
 846         B[i] = gsl_vector_get(x, i*2+1);
 847     }
 848     C = gsl_vector_get(x, nexp*2);
 849
 850     for (i = 0; i < n; i++)
 851     {
 852         t  = i*BD->tDelta;
 853         ct = 0;
 854         for (j = 0; j < nexp; j++)
 855         {
 856             ct += A[j] * exp(B[j] * t);
 857         }
 858         ct += C;
 859         gsl_vector_set (f, i, ct - y[i]);
 860     }
 861     return GSL_SUCCESS;
 862 }
 863
 864 /* The derivative stored in jacobian form (J)*/
 865 static int balFunc_df(const gsl_vector *params, void *data, gsl_matrix *J)
 866 {
 867     balData *BD;
 868     size_t   n, i, j;
 869     double  *y, *A, *B, C, /* There are the parameters. */
 870              t;
 871     int      nexp;
 872
 873     BD   = (balData*)data;
 874     n    = BD->n;
 875     y    = BD->y;
 876     nexp = BD->nexp;
 877
 878     snew(A, nexp);
 879     snew(B, nexp);
 880
 881     for (i = 0; i < nexp; i++)
 882     {
 883         A[i] = gsl_vector_get(params, i*2);
 884         B[i] = gsl_vector_get(params, i*2+1);
 885     }
 886     C = gsl_vector_get(params, nexp*2);
 887     for (i = 0; i < n; i++)
 888     {
 889         t = i*BD->tDelta;
 890         for (j = 0; j < nexp; j++)
 891         {
 892             gsl_matrix_set (J, i, j*2,   exp(B[j]*t));        /* df(t)/dA_j */
 893             gsl_matrix_set (J, i, j*2+1, A[j]*t*exp(B[j]*t)); /* df(t)/dB_j */
 894         }
 895         gsl_matrix_set (J, i, nexp*2, 1);                     /* df(t)/dC */
 896     }
 897     return GSL_SUCCESS;
 898 }
 899
 900 /* Calculation of the function and its derivative */
 901 static int balFunc_fdf(const gsl_vector *params, void *data,
 902                        gsl_vector *f, gsl_matrix *J)
 903 {
 904     balFunc_f(params, data, f);
 905     balFunc_df(params, data, J);
 906     return GSL_SUCCESS;
 907 }
 908 #endif /* HAVE_LIBGSL */
 909
 910 /* Removes the ballistic term from the beginning of the ACF,
 911  * just like in Omer's paper.
 912  */
 913 extern void takeAwayBallistic(double *ct, double *t, int len, real tMax, int nexp, gmx_bool bDerivative)
 914 {
 915
 916     /* Use nonlinear regression with GSL instead.
 917      * Fit with 4 exponentials and one constant term,
 918      * subtract the fatest exponential. */
 919
 920     int      nData, i, status, iter;
 921     balData *BD;
 922     double  *guess,           /* Initial guess. */
 923     *A,                       /* The fitted parameters. (A1, B1, A2, B2,... C) */
 924              a[2],
 925              ddt[2];
 926     gmx_bool sorted;
 927     size_t   n;
 928     size_t   p;
 929
 930     nData = 0;
 931     do
 932     {
 933         nData++;
 934     }
 935     while (t[nData] < tMax+t[0] && nData < len);
 936
 937     p = nexp*2+1;            /* Number of parameters. */
 938
 939 #ifdef HAVE_LIBGSL
 940     const gsl_multifit_fdfsolver_type *T
 941         = gsl_multifit_fdfsolver_lmsder;
 942
 943     gsl_multifit_fdfsolver   *s;           /* The solver itself. */
 944     gsl_multifit_function_fdf fitFunction; /* The function to be fitted. */
 945     gsl_matrix               *covar;       /* Covariance matrix for the parameters.
 946                                             * We'll not use the result, though. */
 947     gsl_vector_view           theParams;
 948
 949     nData = 0;
 950     do
 951     {
 952         nData++;
 953     }
 954     while (t[nData] < tMax+t[0] && nData < len);
 955
 956     guess = NULL;
 957     n     = nData;
 958
 959     snew(guess, p);
 960     snew(A, p);
 961     covar = gsl_matrix_alloc (p, p);
 962
 963     /* Set up an initial gess for the parameters.
 964      * The solver is somewhat sensitive to the initial guess,
 965      * but this worked fine for a TIP3P box with -geminate dd
 966      * EDIT: In fact, this seems like a good starting pont for other watermodels too. */
 967     for (i = 0; i < nexp; i++)
 968     {
 969         guess[i*2]   = 0.1;
 970         guess[i*2+1] = -0.5 + (((double)i)/nexp - 0.5)*0.3;
 971     }
 972     guess[nexp * 2] = 0.01;
 973
 974     theParams = gsl_vector_view_array(guess, p);
 975
 976     snew(BD, 1);
 977     BD->n      = n;
 978     BD->y      = ct;
 979     BD->tDelta = t[1]-t[0];
 980     BD->nexp   = nexp;
 981
 982     fitFunction.f      =  &balFunc_f;
 983     fitFunction.df     =  &balFunc_df;
 984     fitFunction.fdf    =  &balFunc_fdf;
 985     fitFunction.n      =  nData;
 986     fitFunction.p      =  p;
 987     fitFunction.params =  BD;
 988
 989     s = gsl_multifit_fdfsolver_alloc (T, nData, p);
 990     if (s == NULL)
 991     {
 992         gmx_fatal(FARGS, "Could not set up the nonlinear solver.");
 993     }
 994
 995     gsl_multifit_fdfsolver_set(s, &fitFunction, &theParams.vector);
 996
 997     /* \=============================================/ */
 998
 999     iter = 0;
1000     do
1001     {
1002         iter++;
1003         status = gsl_multifit_fdfsolver_iterate (s);
1004
1005         if (status)
1006         {
1007             break;
1008         }
1009         status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
1010     }
1011     while (iter < 5000 && status == GSL_CONTINUE);
1012
1013     if (iter == 5000)
1014     {
1015         fprintf(stderr, "The non-linear fitting did not converge in 5000 steps.\n"
1016                 "Check the quality of the fit!\n");
1017     }
1018     else
1019     {
1020         fprintf(stderr, "Non-linear fitting of ballistic term converged in %d steps.\n\n", (int)iter);
1021     }
1022     for (i = 0; i < nexp; i++)
1023     {
1024         fprintf(stdout, "%c * exp(%c * t) + ", 'A'+(char)i*2, 'B'+(char)i*2);
1025     }
1026
1027     fprintf(stdout, "%c\n", 'A'+(char)nexp*2);
1028     fprintf(stdout, "Here are the actual numbers for A-%c:\n", 'A'+nexp*2);
1029
1030     for (i = 0; i < nexp; i++)
1031     {
1032         A[i*2]   = gsl_vector_get(s->x, i*2);
1033         A[i*2+1] = gsl_vector_get(s->x, i*2+1);
1034         fprintf(stdout, " %g*exp(%g * x) +", A[i*2], A[i*2+1]);
1035     }
1036     A[i*2] = gsl_vector_get(s->x, i*2);        /* The last and constant term */
1037     fprintf(stdout, " %g\n", A[i*2]);
1038
1039     fflush(stdout);
1040
1041     /* Implement some check for parameter quality */
1042     for (i = 0; i < nexp; i++)
1043     {
1044         if (A[i*2] < 0 || A[i*2] > 1)
1045         {
1046             fprintf(stderr, "WARNING: ----------------------------------\n"
1047                     " | A coefficient does not lie within [0,1].\n"
1048                     " | This may or may not be a problem.\n"
1049                     " | Double check the quality of the fit!\n");
1050         }
1051         if (A[i*2+1] > 0)
1052         {
1053             fprintf(stderr, "WARNING: ----------------------------------\n"
1054                     " | One factor in the exponent is positive.\n"
1055                     " | This could be a problem if the coefficient\n"
1056                     " | is large. Double check the quality of the fit!\n");
1057         }
1058     }
1059     if (A[i*2] < 0 || A[i*2] > 1)
1060     {
1061         fprintf(stderr, "WARNING: ----------------------------------\n"
1062                 " | The constant term does not lie within [0,1].\n"
1063                 " | This may or may not be a problem.\n"
1064                 " | Double check the quality of the fit!\n");
1065     }
1066
1067     /* Sort the terms */
1068     sorted = (nexp > 1) ?  FALSE : TRUE;
1069     while (!sorted)
1070     {
1071         sorted = TRUE;
1072         for (i = 0; i < nexp-1; i++)
1073         {
1074             ddt[0] = A[i*2] * A[i*2+1];
1075             ddt[1] = A[i*2+2] * A[i*2+3];
1076
1077             if ((bDerivative && (ddt[0] < 0 && ddt[1] < 0 && ddt[0] > ddt[1])) || /* Compare derivative at t=0... */
1078                 (!bDerivative && (A[i*2+1] > A[i*2+3])))                          /* Or just the coefficient in the exponent */
1079             {
1080                 sorted = FALSE;
1081                 a[0]   = A[i*2];   /* coefficient */
1082                 a[1]   = A[i*2+1]; /* parameter in the exponent */
1083
1084                 A[i*2]   = A[i*2+2];
1085                 A[i*2+1] = A[i*2+3];
1086
1087                 A[i*2+2] = a[0];
1088                 A[i*2+3] = a[1];
1089             }
1090         }
1091     }
1092
1093     /* Subtract the fastest component */
1094     fprintf(stdout, "Fastest component is %g * exp(%g * t)\n"
1095             "Subtracting fastest component from ACF.\n", A[0], A[1]);
1096
1097     for (i = 0; i < len; i++)
1098     {
1099         ct[i] = (ct[i] - A[0] * exp(A[1] * i*BD->tDelta)) / (1-A[0]);
1100     }
1101
1102     sfree(guess);
1103     sfree(A);
1104
1105     gsl_multifit_fdfsolver_free(s);
1106     gsl_matrix_free(covar);
1107     fflush(stdout);
1108
1109 #else
1110     /* We have no gsl. */
1111     fprintf(stderr, "Sorry, can't take away ballistic component without gsl. "
1112             "Recompile using --with-gsl.\n");
1113     return;
1114 #endif /* HAVE_LIBGSL */
1115
1116 }
1117
1118 extern void dumpN(const real *e, const int nn, char *fn)
1119 {
1120     /* For debugging only */
1121     int   i;
1122     FILE *f;
1123     char  standardName[] = "Nt.xvg";
1124     if (fn == NULL)
1125     {
1126         fn = standardName;
1127     }
1128
1129     f = fopen(fn, "w");
1130     fprintf(f,
1131             "@ type XY\n"
1132             "@ xaxis label \"Frame\"\n"
1133             "@ yaxis label \"N\"\n"
1134             "@ s0 line type 3\n");
1135
1136     for (i = 0; i < nn; i++)
1137     {
1138         fprintf(f, "%-10i %-g\n", i, e[i]);
1139     }
1140
1141     fclose(f);
1142 }
1143
1144 static real* d2r(const double *d, const int nn)
1145 {
1146     real *r;
1147     int   i;
1148
1149     snew(r, nn);
1150     for (i = 0; i < nn; i++)
1151     {
1152         r[i] = (real)d[i];
1153     }
1154
1155     return r;
1156 }
1157
1158 static void _patchBad(double *ct, int n, double dy)
1159 {
1160     /* Just do lin. interpolation for now. */
1161     int i;
1162
1163     for (i = 1; i < n; i++)
1164     {
1165         ct[i] = ct[0]+i*dy;
1166     }
1167 }
1168
1169 static void patchBadPart(double *ct, int n)
1170 {
1171     _patchBad(ct, n, (ct[n] - ct[0])/n);
1172 }
1173
1174 static void patchBadTail(double *ct, int n)
1175 {
1176     _patchBad(ct+1, n-1, ct[1]-ct[0]);
1177
1178 }
1179
1180 extern void fixGemACF(double *ct, int len)
1181 {
1182     int      i, j, b, e;
1183     gmx_bool bBad;
1184
1185     /* Let's separate two things:
1186      * - identification of bad parts
1187      * - patching of bad parts.
1188      */
1189
1190     b    = 0; /* Start of a bad stretch */
1191     e    = 0; /* End of a bad stretch */
1192     bBad = FALSE;
1193
1194     /* An acf of binary data must be one at t=0. */
1195     if (abs(ct[0]-1.0) > 1e-6)
1196     {
1197         ct[0] = 1.0;
1198         fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0));
1199     }
1200
1201     for (i = 0; i < len; i++)
1202     {
1203
1204 #ifdef HAS_ISFINITE
1205         if (isfinite(ct[i]))
1206 #elif defined(HAS__ISFINITE)
1207         if (_isfinite(ct[i]))
1208 #else
1209         if (1)
1210 #endif
1211         {
1212             if (!bBad)
1213             {
1214                 /* Still on a good stretch. Proceed.*/
1215                 continue;
1216             }
1217
1218             /* Patch up preceding bad stretch. */
1219             if (i == (len-1))
1220             {
1221                 /* It's the tail */
1222                 if (b <= 1)
1223                 {
1224                     gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting.");
1225                 }
1226                 patchBadTail(&(ct[b-2]), (len-b)+1);
1227             }
1228
1229             e = i;
1230             patchBadPart(&(ct[b-1]), (e-b)+1);
1231             bBad = FALSE;
1232         }
1233         else
1234         {
1235             if (!bBad)
1236             {
1237                 b = i;
1238
1239                 bBad = TRUE;
1240             }
1241         }
1242     }
1243 }