2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
43 #include "gromacs/utility/smalloc.h"
44 #include "gromacs/math/vec.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/utility/gmxomp.h"
50 static void missing_code_message()
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");
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().
61 * The parts menetioned in the previous paragraph were contributed under the BSD license.
65 /* This first part is from complex.c which I recieved from Omer Markowitch.
68 * ------------- from complex.c ------------- */
70 /* Complexation of a paired number (r,i) */
71 static gem_complex
gem_cmplx(double r
, double i
)
79 /* Complexation of a real number, x */
80 static gem_complex
gem_c(double x
)
88 /* Magnitude of a complex number z */
89 static double gem_cx_abs(gem_complex z
)
91 return (sqrt(z
.r
*z
.r
+z
.i
*z
.i
));
94 /* Addition of two complex numbers z1 and z2 */
95 static gem_complex
gem_cxadd(gem_complex z1
, gem_complex z2
)
103 /* Addition of a complex number z1 and a real number r */
104 static gem_complex
gem_cxradd(gem_complex z
, double r
)
112 /* Subtraction of two complex numbers z1 and z2 */
113 static gem_complex
gem_cxsub(gem_complex z1
, gem_complex z2
)
121 /* Multiplication of two complex numbers z1 and z2 */
122 static gem_complex
gem_cxmul(gem_complex z1
, gem_complex z2
)
125 value
.r
= z1
.r
*z2
.r
-z1
.i
*z2
.i
;
126 value
.i
= z1
.r
*z2
.i
+z1
.i
*z2
.r
;
130 /* Square of a complex number z */
131 static gem_complex
gem_cxsq(gem_complex z
)
134 value
.r
= z
.r
*z
.r
-z
.i
*z
.i
;
135 value
.i
= z
.r
*z
.i
*2.;
139 /* multiplication of a complex number z and a real number r */
140 static gem_complex
gem_cxrmul(gem_complex z
, double r
)
148 /* Division of two complex numbers z1 and z2 */
149 static gem_complex
gem_cxdiv(gem_complex z1
, gem_complex z2
)
153 num
= z2
.r
*z2
.r
+z2
.i
*z2
.i
;
156 fprintf(stderr
, "ERROR in gem_cxdiv function\n");
158 value
.r
= (z1
.r
*z2
.r
+z1
.i
*z2
.i
)/num
; value
.i
= (z1
.i
*z2
.r
-z1
.r
*z2
.i
)/num
;
162 /* division of a complex z number by a real number x */
163 static gem_complex
gem_cxrdiv(gem_complex z
, double r
)
171 /* division of a real number r by a complex number x */
172 static gem_complex
gem_rxcdiv(double r
, gem_complex z
)
176 f
= r
/(z
.r
*z
.r
+z
.i
*z
.i
);
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
)
188 value
.r
= exp_z_r
*cos(z
.i
);
189 value
.i
= exp_z_r
*sin(z
.i
);
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
)
199 mag2
= z
.r
*z
.r
+z
.i
*z
.i
;
202 fprintf(stderr
, "ERROR in gem_cxlog func\n");
204 value
.r
= log(sqrt(mag2
));
215 value
.i
= atan2(z
.i
, z
.r
);
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
)
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) */
237 /* Complex power of a complex number z1^z2 */
238 static gem_complex
gem_cxdpow(gem_complex z1
, gem_complex z2
)
241 value
= gem_cxdexp(gem_cxmul(gem_cxlog(z1
), z2
));
245 /* ------------ end of complex.c ------------ */
247 /* This next part was derived from cubic.c, also received from Omer Markovitch.
248 * ------------- from cubic.c ------------- */
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
)
254 double t1
, t2
, two_3
, temp
;
255 gem_complex ctemp
, ct3
;
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
;
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.));
263 ctemp
= gem_rxcdiv(-two_3
*t1
/3., ct3
);
264 ctemp
= gem_cxadd(ctemp
, gem_cxrdiv(ct3
, 3.*two_3
));
266 *gam
= gem_cxadd(gem_cmplx(a
/3., 0.), ctemp
);
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);
275 /* ------------ end of cubic.c ------------ */
277 /* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
278 * ------------- from [cr]error.c ------------- */
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 /************************************************************/
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 */
297 static double gem_erf(double x
)
299 double n
, sum
, temp
, exp2
, xm
, x2
, x4
, x6
, x8
, x10
, x12
;
312 for (n
= 1.; n
<= 2000.; n
+= 1.)
314 temp
*= 2.*x2
/(2.*n
+1.);
316 if (fabs(temp
/sum
) < 1.E
-16)
324 fprintf(stderr
, "In Erf calc - iteration exceeds %lg\n", n
);
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
)
335 sum
*= exp2
; /* now sum is erfc(x) */
338 return x
>= 0.0 ? sum
: -sum
;
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
)
349 ans
= t
* exp(-z
*z
- 1.26551223 +
358 t
*0.17087277)))))))));
360 return x
>= 0.0 ? ans
: 2.0-ans
;
363 /* omega(x)=exp(x^2)erfc(x) */
364 static double gem_omega(double x
)
366 double xm
, ans
, xx
, x4
, x6
, x8
, x10
, x12
;
377 ans
= exp(xx
)*gem_erfc(x
);
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
;
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
)
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
;
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
);
412 for (n
= 1.0, temp
= 0.; n
<= 2000.; n
+= 1.0)
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)
430 fprintf(stderr
, "iteration exceeds %lg\n", n
);
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
));
449 /* ------------ end of [cr]error.c ------------ */
451 /*_ REVERSIBLE GEMINATE RECOMBINATION
453 * Here are the functions for reversible geminate recombination. */
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
)
460 fprintf(stdout
, "Diffusion coefficient is %f A^2/ps\n", D
*1e4
);
461 return (double)(D
*1e4
);
465 static double eq10v2(double theoryCt
[], double time
[], int manytimes
,
466 double ka
, double kd
, t_gemParams
*params
)
468 /* Finding the 3 roots */
470 double kD
, D
, r
, a
, b
, c
, tsqrt
, sumimaginary
;
475 part1
, part2
, part3
, part4
;
480 a
= (1.0 + ka
/kD
) * sqrt(D
)/r
;
484 gem_solve(&alpha
, &beta
, &gamma
, a
, b
, c
);
485 /* Finding the 3 roots */
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) */
493 #pragma omp parallel for \
494 private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
495 reduction(+:sumimaginary) \
498 for (i
= 0; i
< manytimes
; i
++)
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
;
510 sumimaginary
+= c4
.i
* c4
.i
;
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
)
520 return (exp(((double)(i
)) * params
->logQuota
) -1);
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
)
533 /* A few hardcoded things here. For now anyway. */
535 /* p->ka_max = 100; */
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
;
550 /* Parameters used by calcsquare(). Better to calculate them
551 * here than in calcsquare every time it's called. */
553 /* p->logAfterTime = logAfterTime; */
554 tDelta
= (t
[len
-1]-t
[0]) / len
;
557 gmx_fatal(FARGS
, "Time between frames is non-positive!");
564 p
->nExpFit
= nBalExp
;
565 /* p->nLin = logAfterTime / tDelta; */
566 p
->nFitPoints
= nFitPoints
;
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"); */
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 */
580 /* p->logMult = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
581 p
->ballistic
= ballistic
;
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. */
591 static real
* d2r(const double *d
, const int nn
);
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
)
600 int nThreads
, i
, iter
, status
, maxiter
;
601 real size
, d2
, tol
, *dumpdata
;
604 char *dumpstr
, dumpname
[128];
606 missing_code_message();
612 /* Removes the ballistic term from the beginning of the ACF,
613 * just like in Omer's paper.
615 extern void takeAwayBallistic(double gmx_unused
*ct
, double *t
, int len
, real tMax
, int nexp
, gmx_bool gmx_unused bDerivative
)
618 /* Fit with 4 exponentials and one constant term,
619 * subtract the fatest exponential. */
621 int nData
, i
, status
, iter
;
623 double *guess
, /* Initial guess. */
624 *A
, /* The fitted parameters. (A1, B1, A2, B2,... C) */
636 while (t
[nData
] < tMax
+t
[0] && nData
< len
);
638 p
= nexp
*2+1; /* Number of parameters. */
640 missing_code_message();
644 extern void dumpN(const real
*e
, const int nn
, char *fn
)
646 /* For debugging only */
649 char standardName
[] = "Nt.xvg";
658 "@ xaxis label \"Frame\"\n"
659 "@ yaxis label \"N\"\n"
660 "@ s0 line type 3\n");
662 for (i
= 0; i
< nn
; i
++)
664 fprintf(f
, "%-10i %-g\n", i
, e
[i
]);
670 static real
* d2r(const double *d
, const int nn
)
676 for (i
= 0; i
< nn
; i
++)
684 static void _patchBad(double *ct
, int n
, double dy
)
686 /* Just do lin. interpolation for now. */
689 for (i
= 1; i
< n
; i
++)
695 static void patchBadPart(double *ct
, int n
)
697 _patchBad(ct
, n
, (ct
[n
] - ct
[0])/n
);
700 static void patchBadTail(double *ct
, int n
)
702 _patchBad(ct
+1, n
-1, ct
[1]-ct
[0]);
706 extern void fixGemACF(double *ct
, int len
)
711 /* Let's separate two things:
712 * - identification of bad parts
713 * - patching of bad parts.
716 b
= 0; /* Start of a bad stretch */
717 e
= 0; /* End of a bad stretch */
720 /* An acf of binary data must be one at t=0. */
721 if (fabs(ct
[0]-1.0) > 1e-6)
724 fprintf(stderr
, "|ct[0]-1.0| = %1.6f. Setting ct[0] to 1.0.\n", fabs(ct
[0]-1.0));
727 for (i
= 0; i
< len
; i
++)
732 #elif defined(HAS__ISFINITE)
733 if (_isfinite(ct
[i
]))
740 /* Still on a good stretch. Proceed.*/
744 /* Patch up preceding bad stretch. */
750 gmx_fatal(FARGS
, "The ACF is mostly NaN or Inf. Aborting.");
752 patchBadTail(&(ct
[b
-2]), (len
-b
)+1);
756 patchBadPart(&(ct
[b
-1]), (e
-b
)+1);