2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2013, 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.
39 gemNULL, gemNONE, gemDD, gemAD, gemAA, gemA4, gemNR
41 static const char *gemType[] = {NULL, "none", "dd", "ad", "aa", "a4", NULL};
43 /* The first few sections of this file contain functions that were adopted,
44 * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
45 * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
46 * This is also the case with the function eq10v2() in geminate.c.
48 * The parts menetioned in the previous paragraph were contributed under a BSD license.
51 /* This first part is derived from complex.c which I recieved from Omer Markowitch.
54 * ------------- from complex.c ------------- */
57 /* definition of PI */
59 #define PI (acos(-1.0))
62 /* definition of the type complex */
69 /* ------------ end of complex.c ------------ */
71 /* This next part was derived from cerror.c and rerror.c,
72 * also received from Omer Markovitch.
73 * ------------- from [cr]error.c ------------- */
76 #define sPI (sqrt(PI))
79 /* ------------ end of [cr]error.c ------------ */
81 /* ///////////////// REVERSIBLE GEMINATE RECOMBINATION ///////////////////
82 * Here follow routines and structs for reversible geminate recombination.
94 /* Used in the rewritten version of Omer's gem. recomb. analysis */
95 double ka, kd; /* -- || -- results[] */
96 double sigma; /* -- || -- sigma */
97 double D; /* The diffusion coefficient */
98 double kD; /* Replaces kD in analyse_corr_gem3d() */
100 /* The following members are for calcsquare() and takeAwayBallistic() */
101 double tDelta; /* Time between frames */
102 /* double logAfterTime; /\* Time after which we do the lsq calculations on a logarithmic timescale. *\/ */
103 int nFitPoints; /* Number of points to take from the ACF for fitting */
104 double begFit; /* Fit from this time (ps) */
105 double endFit; /* Fit up to this time (ps) */
106 /* double logDelta; */
108 /* To get an equal number of points in the lin and log regime,
109 * we'll use logDelta to calculate where to sample the ACF.
110 * if i and j are indices in the linear and log regime, then:
111 * j = Cexp(A(i+nLin)),
112 * where C = (nLin**2 / len) and A = log(len/nLin) / nLin.
113 * This expands to j = (nLin**2 / len) * exp((i+nLin) * log(len/nLin) / nLin).
114 * We'll save part of our brains and some computing time if we pre-calculate
115 * 1) logDelta = log(len/nLin) / nLin
116 * 2) logPF = nLin**2 / len
117 * and get j = logPF * exp((i+nLin) * logDelta). */
119 /* I will redo this for a fit done entirely in log-log.
121 * nFitPoints' = nFitPoints-1
126 * (j'=len <=> i=nFitPoints')
127 * => A=log(len)/nFitPoints'
128 * => j = exp(i*log(len)/(nFitPoints-1)) -1
130 /* #define GETLOGINDEX(i,params) (params)->logPF * exp(((i)+(params)->nLin) * (params)->logDelta)
133 int nLin; /* Number of timepoints in the linear regime */
134 int len; /* Length of time and ct arrays */
135 int nExpFit; /* Number of exponentials to fit */
136 real ballistic; /* Time before which the ballistic term should be fitted */
137 gmx_bool bDt; /* TRUE => use time derivative at time 0
138 * to find fastest component.
139 * FALSE => use coefficient in exponenetial
140 * to find fastest component. */
145 size_t n; /* Number of data points (lin-log) */
146 double *y; /* The datapoints */
147 double *ctTheory; /* The theoretical ct which will be built by gemFunc_f. */
152 double tDelta; /* time difference between subsequent datapoints */
153 size_t nData; /* real size of the data */
155 double *doubleLogTime;
159 extern void takeAwayBallistic(double *ct, double *t,
161 int nexp, gmx_bool bDerivative);
164 extern t_gemParams *init_gemParams(const double sigma, const double D,
165 const real *t, const int len, const int nFitPoints,
166 const real begFit, const real endFit,
167 const real ballistic, const int nBalExp);
169 /* Fit to geminate recombination model.
170 Returns root mean square error of fit. */
171 extern real fitGemRecomb(double *ct, double *time, double **ctFit,
172 const int nData, t_gemParams *params);
174 extern void dumpN(const real *e, const int nn, char *fn);
176 /* Fix NaN that might appear in the theoretical acf. */
177 extern void fixGemACF(double *ct, int len);