src/gromacs/gmxana/geminate.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "typedefs.h"
  41 #include "gromacs/utility/smalloc.h"
  42 #include "vec.h"
  43 #include "geminate.h"
  44
  45 #include "gromacs/utility/gmxomp.h"
  46
  47 static void missing_code_message()
  48 {
  49     fprintf(stderr, "You have requested code to run that is deprecated.\n");
  50     fprintf(stderr, "Revert to an older GROMACS version or help in porting the code.\n");
  51 }
  52
  53 /* The first few sections of this file contain functions that were adopted,
  54  * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
  55  * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
  56  * This is also the case with the function eq10v2().
  57  *
  58  * The parts menetioned in the previous paragraph were contributed under the BSD license.
  59  */
  60
  61
  62 /* This first part is from complex.c which I recieved from Omer Markowitch.
  63  * - Erik Marklund
  64  *
  65  * ------------- from complex.c ------------- */
  66
  67 /* Complexation of a paired number (r,i)                                     */
  68 static gem_complex gem_cmplx(double r, double i)
  69 {
  70     gem_complex value;
  71     value.r = r;
  72     value.i = i;
  73     return value;
  74 }
  75
  76 /* Complexation of a real number, x */
  77 static gem_complex gem_c(double x)
  78 {
  79     gem_complex value;
  80     value.r = x;
  81     value.i = 0;
  82     return value;
  83 }
  84
  85 /* Magnitude of a complex number z                                           */
  86 static double gem_cx_abs(gem_complex z)
  87 {
  88     return (sqrt(z.r*z.r+z.i*z.i));
  89 }
  90
  91 /* Addition of two complex numbers z1 and z2                                 */
  92 static gem_complex gem_cxadd(gem_complex z1, gem_complex z2)
  93 {
  94     gem_complex value;
  95     value.r = z1.r+z2.r;
  96     value.i = z1.i+z2.i;
  97     return value;
  98 }
  99
 100 /* Addition of a complex number z1 and a real number r */
 101 static gem_complex gem_cxradd(gem_complex z, double r)
 102 {
 103     gem_complex value;
 104     value.r = z.r + r;
 105     value.i = z.i;
 106     return value;
 107 }
 108
 109 /* Subtraction of two complex numbers z1 and z2                              */
 110 static gem_complex gem_cxsub(gem_complex z1, gem_complex z2)
 111 {
 112     gem_complex value;
 113     value.r = z1.r-z2.r;
 114     value.i = z1.i-z2.i;
 115     return value;
 116 }
 117
 118 /* Multiplication of two complex numbers z1 and z2                           */
 119 static gem_complex gem_cxmul(gem_complex z1, gem_complex z2)
 120 {
 121     gem_complex value;
 122     value.r = z1.r*z2.r-z1.i*z2.i;
 123     value.i = z1.r*z2.i+z1.i*z2.r;
 124     return value;
 125 }
 126
 127 /* Square of a complex number z                                              */
 128 static gem_complex gem_cxsq(gem_complex z)
 129 {
 130     gem_complex value;
 131     value.r = z.r*z.r-z.i*z.i;
 132     value.i = z.r*z.i*2.;
 133     return value;
 134 }
 135
 136 /* multiplication of a complex number z and a real number r */
 137 static gem_complex gem_cxrmul(gem_complex z, double r)
 138 {
 139     gem_complex value;
 140     value.r = z.r*r;
 141     value.i = z.i*r;
 142     return value;
 143 }
 144
 145 /* Division of two complex numbers z1 and z2                                 */
 146 static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2)
 147 {
 148     gem_complex value;
 149     double      num;
 150     num = z2.r*z2.r+z2.i*z2.i;
 151     if (num == 0.)
 152     {
 153         fprintf(stderr, "ERROR in gem_cxdiv function\n");
 154     }
 155     value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num;
 156     return value;
 157 }
 158
 159 /* division of a complex z number by a real number x */
 160 static gem_complex gem_cxrdiv(gem_complex z, double r)
 161 {
 162     gem_complex value;
 163     value.r = z.r/r;
 164     value.i = z.i/r;
 165     return value;
 166 }
 167
 168 /* division of a real number r by a complex number x */
 169 static gem_complex gem_rxcdiv(double r, gem_complex z)
 170 {
 171     gem_complex value;
 172     double      f;
 173     f       = r/(z.r*z.r+z.i*z.i);
 174     value.r = f*z.r;
 175     value.i = -f*z.i;
 176     return value;
 177 }
 178
 179 /* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
 180 static gem_complex gem_cxdexp(gem_complex z)
 181 {
 182     gem_complex value;
 183     double      exp_z_r;
 184     exp_z_r = exp(z.r);
 185     value.r = exp_z_r*cos(z.i);
 186     value.i = exp_z_r*sin(z.i);
 187     return value;
 188 }
 189
 190 /* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z),                  */
 191 /*                                  where -PI < Arg(z) < PI                  */
 192 static gem_complex gem_cxlog(gem_complex z)
 193 {
 194     gem_complex value;
 195     double      mag2;
 196     mag2 = z.r*z.r+z.i*z.i;
 197     if (mag2 < 0.)
 198     {
 199         fprintf(stderr, "ERROR in gem_cxlog func\n");
 200     }
 201     value.r = log(sqrt(mag2));
 202     if (z.r == 0.)
 203     {
 204         value.i = PI/2.;
 205         if (z.i < 0.)
 206         {
 207             value.i = -value.i;
 208         }
 209     }
 210     else
 211     {
 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     {
 229         value.r = -value.r;
 230     }
 231     return value;
 232 }
 233
 234 /* Complex power of a complex number z1^z2                                   */
 235 static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2)
 236 {
 237     gem_complex value;
 238     value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2));
 239     return value;
 240 }
 241
 242 /* ------------ end of complex.c ------------ */
 243
 244 /* This next part was derived from cubic.c, also received from Omer Markovitch.
 245  * ------------- from cubic.c ------------- */
 246
 247 /* Solver for a cubic equation: x^3-a*x^2+b*x-c=0                            */
 248 static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam,
 249                       double a, double b, double c)
 250 {
 251     double      t1, t2, two_3, temp;
 252     gem_complex ctemp, ct3;
 253
 254     two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c;
 255     temp  = 4.*t1*t1*t1+t2*t2;
 256
 257     ctemp = gem_cmplx(temp, 0.);   ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp));
 258     ct3   = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.));
 259
 260     ctemp = gem_rxcdiv(-two_3*t1/3., ct3);
 261     ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3));
 262
 263     *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp);
 264
 265     ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.))));
 266     ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam));
 267     ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam);
 268     *al   = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
 269     *be   = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
 270 }
 271
 272 /* ------------ end of cubic.c ------------ */
 273
 274 /* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
 275  * ------------- from [cr]error.c ------------- */
 276
 277 /************************************************************/
 278 /* Real valued error function and related functions         */
 279 /* x, y     : real variables                                */
 280 /* erf(x)   : error function                                */
 281 /* erfc(x)  : complementary error function                  */
 282 /* omega(x) : exp(x*x)*erfc(x)                              */
 283 /* W(x,y)   : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
 284 /************************************************************/
 285
 286 /*---------------------------------------------------------------------------*/
 287 /* Utilzed the series approximation of erf(x)                                */
 288 /* Relative error=|err(x)|/erf(x)<EPS                                        */
 289 /* Handbook of Mathematical functions, Abramowitz, p 297                     */
 290 /* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
 291 /*---------------------------------------------------------------------------*/
 292 /*  This gives erfc(z) correctly only upto >10-15 */
 293
 294 static double gem_erf(double x)
 295 {
 296     double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12;
 297     temp = x;
 298     sum  = temp;
 299     xm   = 26.;
 300     x2   = x*x;
 301     x4   = x2*x2;
 302     x6   = x4*x2;
 303     x8   = x6*x2;
 304     x10  = x8*x2;
 305     x12  = x10*x2;
 306     exp2 = exp(-x2);
 307     if (x <= xm)
 308     {
 309         for (n = 1.; n <= 2000.; n += 1.)
 310         {
 311             temp *= 2.*x2/(2.*n+1.);
 312             sum  += temp;
 313             if (fabs(temp/sum) < 1.E-16)
 314             {
 315                 break;
 316             }
 317         }
 318
 319         if (n >= 2000.)
 320         {
 321             fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n);
 322         }
 323         sum *= 2./sPI*exp2;
 324     }
 325     else
 326     {
 327         /* from the asymptotic expansion of experfc(x) */
 328         sum = (1. - 0.5/x2 + 0.75/x4
 329                - 1.875/x6 + 6.5625/x8
 330                - 29.53125/x10 + 162.421875/x12)
 331             / sPI/x;
 332         sum *= exp2; /* now sum is erfc(x) */
 333         sum  = -sum+1.;
 334     }
 335     return x >= 0.0 ? sum : -sum;
 336 }
 337
 338 /* Result --> Alex's code for erfc and experfc behaves better than this      */
 339 /* complementray error function.                Returns 1.-erf(x)            */
 340 static double gem_erfc(double x)
 341 {
 342     double t, z, ans;
 343     z = fabs(x);
 344     t = 1.0/(1.0+0.5*z);
 345
 346     ans = t * exp(-z*z - 1.26551223 +
 347                   t*(1.00002368 +
 348                      t*(0.37409196 +
 349                         t*(0.09678418 +
 350                            t*(-0.18628806 +
 351                               t*(0.27886807 +
 352                                  t*(-1.13520398 +
 353                                     t*(1.48851587 +
 354                                        t*(-0.82215223 +
 355                                           t*0.17087277)))))))));
 356
 357     return x >= 0.0 ? ans : 2.0-ans;
 358 }
 359
 360 /* omega(x)=exp(x^2)erfc(x)                                                  */
 361 static double gem_omega(double x)
 362 {
 363     double xm, ans, xx, x4, x6, x8, x10, x12;
 364     xm  = 26;
 365     xx  = x*x;
 366     x4  = xx*xx;
 367     x6  = x4*xx;
 368     x8  = x6*xx;
 369     x10 = x8*xx;
 370     x12 = x10*xx;
 371
 372     if (x <= xm)
 373     {
 374         ans = exp(xx)*gem_erfc(x);
 375     }
 376     else
 377     {
 378         /* Asymptotic expansion */
 379         ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x;
 380     }
 381     return ans;
 382 }
 383
 384 /*---------------------------------------------------------------------------*/
 385 /* Utilzed the series approximation of erf(z=x+iy)                           */
 386 /* Relative error=|err(z)|/|erf(z)|<EPS                                      */
 387 /* Handbook of Mathematical functions, Abramowitz, p 299                     */
 388 /* comega(z=x+iy)=cexp(z^2)*cerfc(z)                                         */
 389 /*---------------------------------------------------------------------------*/
 390 static gem_complex gem_comega(gem_complex z)
 391 {
 392     gem_complex value;
 393     double      x, y;
 394     double      sumr, sumi, n, n2, f, temp, temp1;
 395     double      x2, y2, cos_2xy, sin_2xy, cosh_2xy, sinh_2xy, cosh_ny, sinh_ny, exp_y2;
 396
 397     x        = z.r;
 398     y        = z.i;
 399     x2       = x*x;
 400     y2       = y*y;
 401     sumr     = 0.;
 402     sumi     = 0.;
 403     cos_2xy  = cos(2.*x*y);
 404     sin_2xy  = sin(2.*x*y);
 405     cosh_2xy = cosh(2.*x*y);
 406     sinh_2xy = sinh(2.*x*y);
 407     exp_y2   = exp(-y2);
 408
 409     for (n = 1.0, temp = 0.; n <= 2000.; n += 1.0)
 410     {
 411         n2      = n*n;
 412         cosh_ny = cosh(n*y);
 413         sinh_ny = sinh(n*y);
 414         f       = exp(-n2/4.)/(n2+4.*x2);
 415         /* if(f<1.E-200) break; */
 416         sumr += (2.*x - 2.*x*cosh_ny*cos_2xy + n*sinh_ny*sin_2xy)*f;
 417         sumi += (2.*x*cosh_ny*sin_2xy + n*sinh_ny*cos_2xy)*f;
 418         temp1 = sqrt(sumr*sumr+sumi*sumi);
 419         if (fabs((temp1-temp)/temp1) < 1.E-16)
 420         {
 421             break;
 422         }
 423         temp = temp1;
 424     }
 425     if (n == 2000.)
 426     {
 427         fprintf(stderr, "iteration exceeds %lg\n", n);
 428     }
 429     sumr *= 2./PI;
 430     sumi *= 2./PI;
 431
 432     if (x != 0.)
 433     {
 434         f = 1./2./PI/x;
 435     }
 436     else
 437     {
 438         f = 0.;
 439     }
 440     value.r = gem_omega(x)-(f*(1.-cos_2xy)+sumr);
 441     value.i = -(f*sin_2xy+sumi);
 442     value   = gem_cxmul(value, gem_cmplx(exp_y2*cos_2xy, exp_y2*sin_2xy));
 443     return (value);
 444 }
 445
 446 /* ------------ end of [cr]error.c ------------ */
 447
 448 /*_ REVERSIBLE GEMINATE RECOMBINATION
 449  *
 450  * Here are the functions for reversible geminate recombination. */
 451
 452 /* Changes the unit from square cm per s to square Ångström per ps,
 453  * since Omers code uses the latter units while g_mds outputs the former.
 454  * g_hbond expects a diffusion coefficent given in square cm per s. */
 455 static double sqcm_per_s_to_sqA_per_ps (real D)
 456 {
 457     fprintf(stdout, "Diffusion coefficient is %f A^2/ps\n", D*1e4);
 458     return (double)(D*1e4);
 459 }
 460
 461
 462 static double eq10v2(double theoryCt[], double time[], int manytimes,
 463                      double ka, double kd, t_gemParams *params)
 464 {
 465     /* Finding the 3 roots */
 466     int    i;
 467     double kD, D, r, a, b, c, tsqrt, sumimaginary;
 468     gem_complex
 469            alpha, beta, gamma,
 470         c1, c2, c3, c4,
 471         oma, omb, omc,
 472         part1, part2, part3, part4;
 473
 474     kD = params->kD;
 475     D  = params->D;
 476     r  = params->sigma;
 477     a  = (1.0 + ka/kD) * sqrt(D)/r;
 478     b  = kd;
 479     c  = kd * sqrt(D)/r;
 480
 481     gem_solve(&alpha, &beta, &gamma, a, b, c);
 482     /* Finding the 3 roots */
 483
 484     sumimaginary = 0;
 485     part1        = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta,  gamma), gem_cxsub(beta,  gamma)));                 /* 1(2+3)(2-3) */
 486     part2        = gem_cxmul(beta,  gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha)));                 /* 2(3+1)(3-1) */
 487     part3        = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta)));                   /* 3(1+2)(1-2) */
 488     part4        = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */
 489
 490 #pragma omp parallel for \
 491     private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
 492     reduction(+:sumimaginary) \
 493     default(shared) \
 494     schedule(guided)
 495     for (i = 0; i < manytimes; i++)
 496     {
 497         tsqrt         = sqrt(time[i]);
 498         oma           = gem_comega(gem_cxrmul(alpha, tsqrt));
 499         c1            = gem_cxmul(oma, gem_cxdiv(part1, part4));
 500         omb           = gem_comega(gem_cxrmul(beta, tsqrt));
 501         c2            = gem_cxmul(omb, gem_cxdiv(part2, part4));
 502         omc           = gem_comega(gem_cxrmul(gamma, tsqrt));
 503         c3            = gem_cxmul(omc, gem_cxdiv(part3, part4));
 504         c4.r          = c1.r+c2.r+c3.r;
 505         c4.i          = c1.i+c2.i+c3.i;
 506         theoryCt[i]   = c4.r;
 507         sumimaginary += c4.i * c4.i;
 508     }
 509
 510     return sumimaginary;
 511
 512 } /* eq10v2 */
 513
 514 /* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
 515 static double getLogIndex(const int i, const t_gemParams *params)
 516 {
 517     return (exp(((double)(i)) * params->logQuota) -1);
 518 }
 519
 520 extern t_gemParams *init_gemParams(const double sigma, const double D,
 521                                    const real *t, const int len, const int nFitPoints,
 522                                    const real begFit, const real endFit,
 523                                    const real ballistic, const int nBalExp)
 524 {
 525     double       tDelta;
 526     t_gemParams *p;
 527
 528     snew(p, 1);
 529
 530     /* A few hardcoded things here. For now anyway. */
 531 /*   p->ka_min   = 0; */
 532 /*   p->ka_max   = 100; */
 533 /*   p->dka      = 10; */
 534 /*   p->kd_min   = 0; */
 535 /*   p->kd_max   = 2; */
 536 /*   p->dkd      = 0.2; */
 537     p->ka       = 0;
 538     p->kd       = 0;
 539 /*   p->lsq      = -1; */
 540 /*   p->lifetime = 0; */
 541     p->sigma    = sigma*10; /* Omer uses Å, not nm */
 542 /*   p->lsq_old  = 99999; */
 543     p->D        = sqcm_per_s_to_sqA_per_ps(D);
 544     p->kD       = 4*acos(-1.0)*sigma*p->D;
 545
 546
 547     /* Parameters used by calcsquare(). Better to calculate them
 548      * here than in calcsquare every time it's called. */
 549     p->len = len;
 550 /*   p->logAfterTime = logAfterTime; */
 551     tDelta       = (t[len-1]-t[0]) / len;
 552     if (tDelta <= 0)
 553     {
 554         gmx_fatal(FARGS, "Time between frames is non-positive!");
 555     }
 556     else
 557     {
 558         p->tDelta = tDelta;
 559     }
 560
 561     p->nExpFit      = nBalExp;
 562 /*   p->nLin         = logAfterTime / tDelta; */
 563     p->nFitPoints   = nFitPoints;
 564     p->begFit       = begFit;
 565     p->endFit       = endFit;
 566     p->logQuota     = (double)(log(p->len))/(p->nFitPoints-1);
 567 /*   if (p->nLin <= 0) { */
 568 /*     fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
 569 /*     sfree(p); */
 570 /*     return NULL; */
 571 /*   } */
 572 /* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
 573 /* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
 574 /*   p->logPF    = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
 575 /* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */
 576
 577 /*   p->logMult      = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
 578     p->ballistic    =  ballistic;
 579     return p;
 580 }
 581
 582 /* There was a misunderstanding regarding the fitting. From our
 583  * recent correspondence it appears that Omer's code require
 584  * the ACF data on a log-scale and does not operate on the raw data.
 585  * This needs to be redone in gemFunc_residual() as well as in the
 586  * t_gemParams structure. */
 587
 588 static real* d2r(const double *d, const int nn);
 589
 590 extern real fitGemRecomb(double gmx_unused      *ct,
 591                          double gmx_unused      *time,
 592                          double gmx_unused     **ctFit,
 593                          const int gmx_unused    nData,
 594                          t_gemParams gmx_unused *params)
 595 {
 596
 597     int         nThreads, i, iter, status, maxiter;
 598     real        size, d2, tol, *dumpdata;
 599     size_t      p, n;
 600     gemFitData *GD;
 601     char       *dumpstr, dumpname[128];
 602
 603     missing_code_message();
 604     return -1;
 605
 606 }
 607
 608
 609 /* Removes the ballistic term from the beginning of the ACF,
 610  * just like in Omer's paper.
 611  */
 612 extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative)
 613 {
 614
 615     /* Fit with 4 exponentials and one constant term,
 616      * subtract the fatest exponential. */
 617
 618     int      nData, i, status, iter;
 619     balData *BD;
 620     double  *guess,           /* Initial guess. */
 621     *A,                       /* The fitted parameters. (A1, B1, A2, B2,... C) */
 622              a[2],
 623              ddt[2];
 624     gmx_bool sorted;
 625     size_t   n;
 626     size_t   p;
 627
 628     nData = 0;
 629     do
 630     {
 631         nData++;
 632     }
 633     while (t[nData] < tMax+t[0] && nData < len);
 634
 635     p = nexp*2+1;            /* Number of parameters. */
 636
 637     missing_code_message();
 638     return;
 639 }
 640
 641 extern void dumpN(const real *e, const int nn, char *fn)
 642 {
 643     /* For debugging only */
 644     int   i;
 645     FILE *f;
 646     char  standardName[] = "Nt.xvg";
 647     if (fn == NULL)
 648     {
 649         fn = standardName;
 650     }
 651
 652     f = fopen(fn, "w");
 653     fprintf(f,
 654             "@ type XY\n"
 655             "@ xaxis label \"Frame\"\n"
 656             "@ yaxis label \"N\"\n"
 657             "@ s0 line type 3\n");
 658
 659     for (i = 0; i < nn; i++)
 660     {
 661         fprintf(f, "%-10i %-g\n", i, e[i]);
 662     }
 663
 664     fclose(f);
 665 }
 666
 667 static real* d2r(const double *d, const int nn)
 668 {
 669     real *r;
 670     int   i;
 671
 672     snew(r, nn);
 673     for (i = 0; i < nn; i++)
 674     {
 675         r[i] = (real)d[i];
 676     }
 677
 678     return r;
 679 }
 680
 681 static void _patchBad(double *ct, int n, double dy)
 682 {
 683     /* Just do lin. interpolation for now. */
 684     int i;
 685
 686     for (i = 1; i < n; i++)
 687     {
 688         ct[i] = ct[0]+i*dy;
 689     }
 690 }
 691
 692 static void patchBadPart(double *ct, int n)
 693 {
 694     _patchBad(ct, n, (ct[n] - ct[0])/n);
 695 }
 696
 697 static void patchBadTail(double *ct, int n)
 698 {
 699     _patchBad(ct+1, n-1, ct[1]-ct[0]);
 700
 701 }
 702
 703 extern void fixGemACF(double *ct, int len)
 704 {
 705     int      i, j, b, e;
 706     gmx_bool bBad;
 707
 708     /* Let's separate two things:
 709      * - identification of bad parts
 710      * - patching of bad parts.
 711      */
 712
 713     b    = 0; /* Start of a bad stretch */
 714     e    = 0; /* End of a bad stretch */
 715     bBad = FALSE;
 716
 717     /* An acf of binary data must be one at t=0. */
 718     if (abs(ct[0]-1.0) > 1e-6)
 719     {
 720         ct[0] = 1.0;
 721         fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0));
 722     }
 723
 724     for (i = 0; i < len; i++)
 725     {
 726
 727 #ifdef HAS_ISFINITE
 728         if (isfinite(ct[i]))
 729 #elif defined(HAS__ISFINITE)
 730         if (_isfinite(ct[i]))
 731 #else
 732         if (1)
 733 #endif
 734         {
 735             if (!bBad)
 736             {
 737                 /* Still on a good stretch. Proceed.*/
 738                 continue;
 739             }
 740
 741             /* Patch up preceding bad stretch. */
 742             if (i == (len-1))
 743             {
 744                 /* It's the tail */
 745                 if (b <= 1)
 746                 {
 747                     gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting.");
 748                 }
 749                 patchBadTail(&(ct[b-2]), (len-b)+1);
 750             }
 751
 752             e = i;
 753             patchBadPart(&(ct[b-1]), (e-b)+1);
 754             bBad = FALSE;
 755         }
 756         else
 757         {
 758             if (!bBad)
 759             {
 760                 b = i;
 761
 762                 bBad = TRUE;
 763             }
 764         }
 765     }
 766 }