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