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