X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=blobdiff_plain;f=src%2Fgromacs%2Fgmxana%2Fgeminate.c;fp=src%2Fgromacs%2Fgmxana%2Fgeminate.c;h=0000000000000000000000000000000000000000;hb=aef61d9804d107c75e373f3e4c29db5c8f18d6a4;hp=d06a4b41d8e91e8ffde3df410454fe98a4181083;hpb=1f61b284053fa52ce8d076fcefb2bffeafd50042;p=alexxy%2Fgromacs.git diff --git a/src/gromacs/gmxana/geminate.c b/src/gromacs/gmxana/geminate.c deleted file mode 100644 index d06a4b41d8..0000000000 --- a/src/gromacs/gmxana/geminate.c +++ /dev/null @@ -1,767 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by - * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, - * and including many others, as listed in the AUTHORS file in the - * top-level source directory and at http://www.gromacs.org. - * - * GROMACS is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public License - * as published by the Free Software Foundation; either version 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -#include "gmxpre.h" - -#include "geminate.h" - -#include -#include - -#include "gromacs/legacyheaders/typedefs.h" -#include "gromacs/math/vec.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/gmxomp.h" -#include "gromacs/utility/smalloc.h" - -static void missing_code_message() -{ - fprintf(stderr, "You have requested code to run that is deprecated.\n"); - fprintf(stderr, "Revert to an older GROMACS version or help in porting the code.\n"); -} - -/* The first few sections of this file contain functions that were adopted, - * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se, - * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url). - * This is also the case with the function eq10v2(). - * - * The parts menetioned in the previous paragraph were contributed under the BSD license. - */ - - -/* This first part is from complex.c which I recieved from Omer Markowitch. - * - Erik Marklund - * - * ------------- from complex.c ------------- */ - -/* Complexation of a paired number (r,i) */ -static gem_complex gem_cmplx(double r, double i) -{ - gem_complex value; - value.r = r; - value.i = i; - return value; -} - -/* Complexation of a real number, x */ -static gem_complex gem_c(double x) -{ - gem_complex value; - value.r = x; - value.i = 0; - return value; -} - -/* Magnitude of a complex number z */ -static double gem_cx_abs(gem_complex z) -{ - return (sqrt(z.r*z.r+z.i*z.i)); -} - -/* Addition of two complex numbers z1 and z2 */ -static gem_complex gem_cxadd(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r+z2.r; - value.i = z1.i+z2.i; - return value; -} - -/* Addition of a complex number z1 and a real number r */ -static gem_complex gem_cxradd(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r + r; - value.i = z.i; - return value; -} - -/* Subtraction of two complex numbers z1 and z2 */ -static gem_complex gem_cxsub(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r-z2.r; - value.i = z1.i-z2.i; - return value; -} - -/* Multiplication of two complex numbers z1 and z2 */ -static gem_complex gem_cxmul(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value.r = z1.r*z2.r-z1.i*z2.i; - value.i = z1.r*z2.i+z1.i*z2.r; - return value; -} - -/* Square of a complex number z */ -static gem_complex gem_cxsq(gem_complex z) -{ - gem_complex value; - value.r = z.r*z.r-z.i*z.i; - value.i = z.r*z.i*2.; - return value; -} - -/* multiplication of a complex number z and a real number r */ -static gem_complex gem_cxrmul(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r*r; - value.i = z.i*r; - return value; -} - -/* Division of two complex numbers z1 and z2 */ -static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2) -{ - gem_complex value; - double num; - num = z2.r*z2.r+z2.i*z2.i; - if (num == 0.) - { - fprintf(stderr, "ERROR in gem_cxdiv function\n"); - } - value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num; - return value; -} - -/* division of a complex z number by a real number x */ -static gem_complex gem_cxrdiv(gem_complex z, double r) -{ - gem_complex value; - value.r = z.r/r; - value.i = z.i/r; - return value; -} - -/* division of a real number r by a complex number x */ -static gem_complex gem_rxcdiv(double r, gem_complex z) -{ - gem_complex value; - double f; - f = r/(z.r*z.r+z.i*z.i); - value.r = f*z.r; - value.i = -f*z.i; - return value; -} - -/* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/ -static gem_complex gem_cxdexp(gem_complex z) -{ - gem_complex value; - double exp_z_r; - exp_z_r = exp(z.r); - value.r = exp_z_r*cos(z.i); - value.i = exp_z_r*sin(z.i); - return value; -} - -/* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z), */ -/* where -PI < Arg(z) < PI */ -static gem_complex gem_cxlog(gem_complex z) -{ - gem_complex value; - double mag2; - mag2 = z.r*z.r+z.i*z.i; - if (mag2 < 0.) - { - fprintf(stderr, "ERROR in gem_cxlog func\n"); - } - value.r = log(sqrt(mag2)); - if (z.r == 0.) - { - value.i = PI/2.; - if (z.i < 0.) - { - value.i = -value.i; - } - } - else - { - value.i = atan2(z.i, z.r); - } - return value; -} - -/* Square root of a complex number z = |z| exp(I*the) -- z^(1/2) */ -/* z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */ -/* where 0 < the < 2*PI */ -static gem_complex gem_cxdsqrt(gem_complex z) -{ - gem_complex value; - double sq; - sq = gem_cx_abs(z); - value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */ - value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */ - if (z.i < 0.) - { - value.r = -value.r; - } - return value; -} - -/* Complex power of a complex number z1^z2 */ -static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2) -{ - gem_complex value; - value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2)); - return value; -} - -/* ------------ end of complex.c ------------ */ - -/* This next part was derived from cubic.c, also received from Omer Markovitch. - * ------------- from cubic.c ------------- */ - -/* Solver for a cubic equation: x^3-a*x^2+b*x-c=0 */ -static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam, - double a, double b, double c) -{ - double t1, t2, two_3, temp; - gem_complex ctemp, ct3; - - two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c; - temp = 4.*t1*t1*t1+t2*t2; - - ctemp = gem_cmplx(temp, 0.); ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp)); - ct3 = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.)); - - ctemp = gem_rxcdiv(-two_3*t1/3., ct3); - ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3)); - - *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp); - - ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.)))); - ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam)); - ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam); - *al = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5); - *be = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5); -} - -/* ------------ end of cubic.c ------------ */ - -/* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch. - * ------------- from [cr]error.c ------------- */ - -/************************************************************/ -/* Real valued error function and related functions */ -/* x, y : real variables */ -/* erf(x) : error function */ -/* erfc(x) : complementary error function */ -/* omega(x) : exp(x*x)*erfc(x) */ -/* W(x,y) : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */ -/************************************************************/ - -/*---------------------------------------------------------------------------*/ -/* Utilzed the series approximation of erf(x) */ -/* Relative error=|err(x)|/erf(x)=6 series sum deteriorates -> Asymptotic series used instead */ -/*---------------------------------------------------------------------------*/ -/* This gives erfc(z) correctly only upto >10-15 */ - -static double gem_erf(double x) -{ - double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12; - temp = x; - sum = temp; - xm = 26.; - x2 = x*x; - x4 = x2*x2; - x6 = x4*x2; - x8 = x6*x2; - x10 = x8*x2; - x12 = x10*x2; - exp2 = exp(-x2); - if (x <= xm) - { - for (n = 1.; n <= 2000.; n += 1.) - { - temp *= 2.*x2/(2.*n+1.); - sum += temp; - if (fabs(temp/sum) < 1.E-16) - { - break; - } - } - - if (n >= 2000.) - { - fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n); - } - sum *= 2./sPI*exp2; - } - else - { - /* from the asymptotic expansion of experfc(x) */ - sum = (1. - 0.5/x2 + 0.75/x4 - - 1.875/x6 + 6.5625/x8 - - 29.53125/x10 + 162.421875/x12) - / sPI/x; - sum *= exp2; /* now sum is erfc(x) */ - sum = -sum+1.; - } - return x >= 0.0 ? sum : -sum; -} - -/* Result --> Alex's code for erfc and experfc behaves better than this */ -/* complementray error function. Returns 1.-erf(x) */ -static double gem_erfc(double x) -{ - double t, z, ans; - z = fabs(x); - t = 1.0/(1.0+0.5*z); - - ans = t * exp(-z*z - 1.26551223 + - t*(1.00002368 + - t*(0.37409196 + - t*(0.09678418 + - t*(-0.18628806 + - t*(0.27886807 + - t*(-1.13520398 + - t*(1.48851587 + - t*(-0.82215223 + - t*0.17087277))))))))); - - return x >= 0.0 ? ans : 2.0-ans; -} - -/* omega(x)=exp(x^2)erfc(x) */ -static double gem_omega(double x) -{ - double xm, ans, xx, x4, x6, x8, x10, x12; - xm = 26; - xx = x*x; - x4 = xx*xx; - x6 = x4*xx; - x8 = x6*xx; - x10 = x8*xx; - x12 = x10*xx; - - if (x <= xm) - { - ans = exp(xx)*gem_erfc(x); - } - else - { - /* Asymptotic expansion */ - ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x; - } - return ans; -} - -/*---------------------------------------------------------------------------*/ -/* Utilzed the series approximation of erf(z=x+iy) */ -/* Relative error=|err(z)|/|erf(z)|kD; - D = params->D; - r = params->sigma; - a = (1.0 + ka/kD) * sqrt(D)/r; - b = kd; - c = kd * sqrt(D)/r; - - gem_solve(&alpha, &beta, &gamma, a, b, c); - /* Finding the 3 roots */ - - sumimaginary = 0; - part1 = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta, gamma), gem_cxsub(beta, gamma))); /* 1(2+3)(2-3) */ - part2 = gem_cxmul(beta, gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha))); /* 2(3+1)(3-1) */ - part3 = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta))); /* 3(1+2)(1-2) */ - part4 = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */ - -#pragma omp parallel for \ - private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \ - reduction(+:sumimaginary) \ - default(shared) \ - schedule(guided) - for (i = 0; i < manytimes; i++) - { - tsqrt = sqrt(time[i]); - oma = gem_comega(gem_cxrmul(alpha, tsqrt)); - c1 = gem_cxmul(oma, gem_cxdiv(part1, part4)); - omb = gem_comega(gem_cxrmul(beta, tsqrt)); - c2 = gem_cxmul(omb, gem_cxdiv(part2, part4)); - omc = gem_comega(gem_cxrmul(gamma, tsqrt)); - c3 = gem_cxmul(omc, gem_cxdiv(part3, part4)); - c4.r = c1.r+c2.r+c3.r; - c4.i = c1.i+c2.i+c3.i; - theoryCt[i] = c4.r; - sumimaginary += c4.i * c4.i; - } - - return sumimaginary; - -} /* eq10v2 */ - -/* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */ -static double getLogIndex(const int i, const t_gemParams *params) -{ - return gmx_expm1(((double)(i)) * params->logQuota); -} - -extern t_gemParams *init_gemParams(const double sigma, const double D, - const real *t, const int len, const int nFitPoints, - const real begFit, const real endFit, - const real ballistic, const int nBalExp) -{ - double tDelta; - t_gemParams *p; - - snew(p, 1); - - /* A few hardcoded things here. For now anyway. */ -/* p->ka_min = 0; */ -/* p->ka_max = 100; */ -/* p->dka = 10; */ -/* p->kd_min = 0; */ -/* p->kd_max = 2; */ -/* p->dkd = 0.2; */ - p->ka = 0; - p->kd = 0; -/* p->lsq = -1; */ -/* p->lifetime = 0; */ - p->sigma = sigma*10; /* Omer uses Å, not nm */ -/* p->lsq_old = 99999; */ - p->D = sqcm_per_s_to_sqA_per_ps(D); - p->kD = 4*acos(-1.0)*sigma*p->D; - - - /* Parameters used by calcsquare(). Better to calculate them - * here than in calcsquare every time it's called. */ - p->len = len; -/* p->logAfterTime = logAfterTime; */ - tDelta = (t[len-1]-t[0]) / len; - if (tDelta <= 0) - { - gmx_fatal(FARGS, "Time between frames is non-positive!"); - } - else - { - p->tDelta = tDelta; - } - - p->nExpFit = nBalExp; -/* p->nLin = logAfterTime / tDelta; */ - p->nFitPoints = nFitPoints; - p->begFit = begFit; - p->endFit = endFit; - p->logQuota = (double)(log(p->len))/(p->nFitPoints-1); -/* if (p->nLin <= 0) { */ -/* fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */ -/* sfree(p); */ -/* return NULL; */ -/* } */ -/* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */ -/* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */ -/* p->logPF = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */ -/* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */ - -/* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */ - p->ballistic = ballistic; - return p; -} - -/* There was a misunderstanding regarding the fitting. From our - * recent correspondence it appears that Omer's code require - * the ACF data on a log-scale and does not operate on the raw data. - * This needs to be redone in gemFunc_residual() as well as in the - * t_gemParams structure. */ - -static real* d2r(const double *d, const int nn); - -extern real fitGemRecomb(double gmx_unused *ct, - double gmx_unused *time, - double gmx_unused **ctFit, - const int gmx_unused nData, - t_gemParams gmx_unused *params) -{ - - int nThreads, i, iter, status, maxiter; - real size, d2, tol, *dumpdata; - size_t p, n; - gemFitData *GD; - char *dumpstr, dumpname[128]; - - missing_code_message(); - return -1; - -} - - -/* Removes the ballistic term from the beginning of the ACF, - * just like in Omer's paper. - */ -extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative) -{ - - /* Fit with 4 exponentials and one constant term, - * subtract the fatest exponential. */ - - int nData, i, status, iter; - balData *BD; - double *guess, /* Initial guess. */ - *A, /* The fitted parameters. (A1, B1, A2, B2,... C) */ - a[2], - ddt[2]; - gmx_bool sorted; - size_t n; - size_t p; - - nData = 0; - do - { - nData++; - } - while (t[nData] < tMax+t[0] && nData < len); - - p = nexp*2+1; /* Number of parameters. */ - - missing_code_message(); - return; -} - -extern void dumpN(const real *e, const int nn, char *fn) -{ - /* For debugging only */ - int i; - FILE *f; - char standardName[] = "Nt.xvg"; - if (fn == NULL) - { - fn = standardName; - } - - f = fopen(fn, "w"); - fprintf(f, - "@ type XY\n" - "@ xaxis label \"Frame\"\n" - "@ yaxis label \"N\"\n" - "@ s0 line type 3\n"); - - for (i = 0; i < nn; i++) - { - fprintf(f, "%-10i %-g\n", i, e[i]); - } - - fclose(f); -} - -static real* d2r(const double *d, const int nn) -{ - real *r; - int i; - - snew(r, nn); - for (i = 0; i < nn; i++) - { - r[i] = (real)d[i]; - } - - return r; -} - -static void _patchBad(double *ct, int n, double dy) -{ - /* Just do lin. interpolation for now. */ - int i; - - for (i = 1; i < n; i++) - { - ct[i] = ct[0]+i*dy; - } -} - -static void patchBadPart(double *ct, int n) -{ - _patchBad(ct, n, (ct[n] - ct[0])/n); -} - -static void patchBadTail(double *ct, int n) -{ - _patchBad(ct+1, n-1, ct[1]-ct[0]); - -} - -extern void fixGemACF(double *ct, int len) -{ - int i, j, b, e; - gmx_bool bBad; - - /* Let's separate two things: - * - identification of bad parts - * - patching of bad parts. - */ - - b = 0; /* Start of a bad stretch */ - e = 0; /* End of a bad stretch */ - bBad = FALSE; - - /* An acf of binary data must be one at t=0. */ - if (fabs(ct[0]-1.0) > 1e-6) - { - ct[0] = 1.0; - fprintf(stderr, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct[0]-1.0)); - } - - for (i = 0; i < len; i++) - { - -#ifdef HAS_ISFINITE - if (isfinite(ct[i])) -#elif defined(HAS__ISFINITE) - if (_isfinite(ct[i])) -#else - if (1) -#endif - { - if (!bBad) - { - /* Still on a good stretch. Proceed.*/ - continue; - } - - /* Patch up preceding bad stretch. */ - if (i == (len-1)) - { - /* It's the tail */ - if (b <= 1) - { - gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting."); - } - patchBadTail(&(ct[b-2]), (len-b)+1); - } - - e = i; - patchBadPart(&(ct[b-1]), (e-b)+1); - bBad = FALSE; - } - else - { - if (!bBad) - { - b = i; - - bBad = TRUE; - } - } - } -}