| File: | gromacs/gmxana/geminate.c |
| Location: | line 409, column 5 |
| Description: | Value stored to 'sinh_2xy' is never read |
| 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_H1 |
| 36 | #include <config.h> |
| 37 | #endif |
| 38 | |
| 39 | #include <math.h> |
| 40 | #include <stdlib.h> |
| 41 | |
| 42 | #include "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(stderrstderr, "You have requested code to run that is deprecated.\n"); |
| 53 | fprintf(stderrstderr, "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(stderrstderr, "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(stderrstderr, "ERROR in gem_cxlog func\n"); |
| 203 | } |
| 204 | value.r = log(sqrt(mag2)); |
| 205 | if (z.r == 0.) |
| 206 | { |
| 207 | value.i = PI(acos(-1.0))/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(stderrstderr, "In Erf calc - iteration exceeds %lg\n", n); |
| 325 | } |
| 326 | sum *= 2./sPI(sqrt((acos(-1.0))))*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(sqrt((acos(-1.0))))/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(sqrt((acos(-1.0))))/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); |
Value stored to 'sinh_2xy' is never read | |
| 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(stderrstderr, "iteration exceeds %lg\n", n); |
| 431 | } |
| 432 | sumr *= 2./PI(acos(-1.0)); |
| 433 | sumi *= 2./PI(acos(-1.0)); |
| 434 | |
| 435 | if (x != 0.) |
| 436 | { |
| 437 | f = 1./2./PI(acos(-1.0))/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(stdoutstdout, "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)(p) = save_calloc("p", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/geminate.c" , 531, (1), sizeof(*(p))); |
| 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(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/geminate.c" , 557, "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__attribute__ ((unused)) *ct, |
| 594 | double gmx_unused__attribute__ ((unused)) *time, |
| 595 | double gmx_unused__attribute__ ((unused)) **ctFit, |
| 596 | const int gmx_unused__attribute__ ((unused)) nData, |
| 597 | t_gemParams gmx_unused__attribute__ ((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__attribute__ ((unused)) *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused__attribute__ ((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((void*)0)) |
| 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)(r) = save_calloc("r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/geminate.c" , 675, (nn), sizeof(*(r))); |
| 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 = FALSE0; |
| 719 | |
| 720 | /* An acf of binary data must be one at t=0. */ |
| 721 | if (abs(ct[0]-1.0) > 1e-6) |
| 722 | { |
| 723 | ct[0] = 1.0; |
| 724 | fprintf(stderrstderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0)); |
| 725 | } |
| 726 | |
| 727 | for (i = 0; i < len; i++) |
| 728 | { |
| 729 | |
| 730 | #ifdef HAS_ISFINITE |
| 731 | if (isfinite(ct[i])(sizeof (ct[i]) == sizeof (float) ? __finitef (ct[i]) : sizeof (ct[i]) == sizeof (double) ? __finite (ct[i]) : __finitel (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(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/geminate.c" , 750, "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 = FALSE0; |
| 758 | } |
| 759 | else |
| 760 | { |
| 761 | if (!bBad) |
| 762 | { |
| 763 | b = i; |
| 764 | |
| 765 | bBad = TRUE1; |
| 766 | } |
| 767 | } |
| 768 | } |
| 769 | } |