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