*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013, by the GROMACS development team, led by
+ * Copyright (c) 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.
* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <math.h>
+#include "gmxpre.h"
-/*#define HAVE_NN_LOOPS*/
+#include "config.h"
+
+#include <math.h>
#include "gromacs/commandline/pargs.h"
-#include "copyrite.h"
-#include "sysstuff.h"
-#include "txtdump.h"
-#include "physics.h"
-#include "macros.h"
-#include "gmx_fatal.h"
-#include "index.h"
-#include "smalloc.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "gstat.h"
-#include "string2.h"
-#include "pbc.h"
-#include "correl.h"
-#include "gmx_ana.h"
-#include "geminate.h"
-
-#include "gromacs/fileio/futil.h"
#include "gromacs/fileio/matio.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/correl.h"
+#include "gromacs/gmxana/geminate.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/gmxana/gstat.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/viewit.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/smalloc.h"
+
+/*#define HAVE_NN_LOOPS*/
typedef short int t_E;
typedef int t_EEst;
{
int i;
- if (ISDON(datable[id]) || !datable)
+ if (!datable || ISDON(datable[id]))
{
if (ddd->dptr[id] == NOTSET) /* New donor */
{
* This could be implemented slightly more efficient, but the code
* would get much more complicated.
*/
-static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+static gmx_inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
{
return ((n == 1) ? x : bTric && bEdge ? 0 : (x-1));
}
-static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+static gmx_inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
{
return ((n == 1) ? x : bTric && bEdge ? (n-1) : (x+1));
}
-static inline int grid_mod(int j, int n)
+static gmx_inline int grid_mod(int j, int n)
{
return (j+n) % (n);
}
integral += x1;
}
integral *= dt;
- ffclose(fp);
+ gmx_ffclose(fp);
printf("%s lifetime = %.2f ps\n", bContact ? "Contact" : "HB", integral);
printf("Note that the lifetime obtained in this manner is close to useless\n");
printf("Use the -ac option instead and check the Forward lifetime\n");
{
return;
}
- fp = ffopen("debug-ac.xvg", "w");
+ fp = gmx_ffopen("debug-ac.xvg", "w");
for (j = 0; (j < nframes); j++)
{
fprintf(fp, "%10.3f", hb->time[j]);
}
fprintf(fp, "\n");
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
static real calc_dg(real tau, real temp)
real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt;
} t_luzar;
-#ifdef HAVE_LIBGSL
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_sf.h>
-#include <gsl/gsl_version.h>
-
-static double my_f(const gsl_vector *v, void *params)
-{
- t_luzar *tl = (t_luzar *)params;
- int i;
- double tol = 1e-16, chi2 = 0;
- double di;
- real k, kp;
-
- for (i = 0; (i < tl->nparams); i++)
- {
- tl->kkk[i] = gsl_vector_get(v, i);
- }
- k = tl->kkk[0];
- kp = tl->kkk[1];
-
- for (i = tl->n0; (i < tl->n1); i += tl->ndelta)
- {
- di = sqr(k*tl->sigma_ct[i]) + sqr(kp*tl->sigma_nt[i]) + sqr(tl->sigma_kt[i]);
- /*di = 1;*/
- if (di > tol)
- {
- chi2 += sqr(k*tl->ct[i]-kp*tl->nt[i]-tl->kt[i])/di;
- }
-
- else
- {
- fprintf(stderr, "WARNING: sigma_ct = %g, sigma_nt = %g, sigma_kt = %g\n"
- "di = %g k = %g kp = %g\n", tl->sigma_ct[i],
- tl->sigma_nt[i], tl->sigma_kt[i], di, k, kp);
- }
- }
-#ifdef DEBUG
- chi2 = 0.3*sqr(k-0.6)+0.7*sqr(kp-1.3);
-#endif
- return chi2;
-}
-
-static real optimize_luzar_parameters(FILE *fp, t_luzar *tl, int maxiter,
- real tol)
-{
- real size, d2;
- int iter = 0;
- int status = 0;
- int i;
-
- const gsl_multimin_fminimizer_type *T;
- gsl_multimin_fminimizer *s;
-
- gsl_vector *x, *dx;
- gsl_multimin_function my_func;
-
- my_func.f = &my_f;
- my_func.n = tl->nparams;
- my_func.params = (void *) tl;
-
- /* Starting point */
- x = gsl_vector_alloc (my_func.n);
- for (i = 0; (i < my_func.n); i++)
- {
- gsl_vector_set (x, i, tl->kkk[i]);
- }
-
- /* Step size, different for each of the parameters */
- dx = gsl_vector_alloc (my_func.n);
- for (i = 0; (i < my_func.n); i++)
- {
- gsl_vector_set (dx, i, 0.01*tl->kkk[i]);
- }
-
- T = gsl_multimin_fminimizer_nmsimplex;
- s = gsl_multimin_fminimizer_alloc (T, my_func.n);
-
- gsl_multimin_fminimizer_set (s, &my_func, x, dx);
- gsl_vector_free (x);
- gsl_vector_free (dx);
-
- if (fp)
- {
- fprintf(fp, "%5s %12s %12s %12s %12s\n", "Iter", "k", "kp", "NM Size", "Chi2");
- }
-
- do
- {
- iter++;
- status = gsl_multimin_fminimizer_iterate (s);
-
- if (status != 0)
- {
- gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
- gsl_multimin_fminimizer_name(s));
- }
-
- d2 = gsl_multimin_fminimizer_minimum(s);
- size = gsl_multimin_fminimizer_size(s);
- status = gsl_multimin_test_size(size, tol);
-
- if (status == GSL_SUCCESS)
- {
- if (fp)
- {
- fprintf(fp, "Minimum found using %s at:\n",
- gsl_multimin_fminimizer_name(s));
- }
- }
-
- if (fp)
- {
- fprintf(fp, "%5d", iter);
- for (i = 0; (i < my_func.n); i++)
- {
- fprintf(fp, " %12.4e", gsl_vector_get (s->x, i));
- }
- fprintf (fp, " %12.4e %12.4e\n", size, d2);
- }
- }
- while ((status == GSL_CONTINUE) && (iter < maxiter));
-
- gsl_multimin_fminimizer_free (s);
-
- return d2;
-}
-
-static real quality_of_fit(real chi2, int N)
-{
- return gsl_sf_gamma_inc_Q((N-2)/2.0, chi2/2.0);
-}
-
-#else
-static real optimize_luzar_parameters(FILE gmx_unused *fp,
- t_luzar gmx_unused *tl,
- int gmx_unused maxiter,
- real gmx_unused tol)
-{
- fprintf(stderr, "This program needs the GNU scientific library to work.\n");
-
- return -1;
-}
-static real quality_of_fit(real gmx_unused chi2, int gmx_unused N)
-{
- fprintf(stderr, "This program needs the GNU scientific library to work.\n");
-
- return -1;
-}
-
-#endif
-
static real compute_weighted_rates(int n, real t[], real ct[], real nt[],
real kt[], real sigma_ct[], real sigma_nt[],
real sigma_kt[], real *k, real *kp,
tl.kkk[0] = *k;
tl.kkk[1] = *kp;
- chi2 = optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
+ chi2 = 0; /*optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
*k = tl.kkk[0];
*kp = tl.kkk[1] = *kp;
tl.ndelta = NK;
for (j = 0; (j < NK); j++)
{
- (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3);
+ /* (void) optimize_luzar_parameters(debug, &tl, 1000, 1e-3); */
kkk += tl.kkk[0];
kkp += tl.kkk[1];
kk2 += sqr(tl.kkk[0]);
chi22 = compute_weighted_rates(n, t, ct, nt, kt, sigma_ct, sigma_nt,
sigma_kt, &k, &kp,
&sigma_k, &sigma_kp, fit_start);
- Q = quality_of_fit(chi2, 2);
+ Q = 0; /* quality_of_fit(chi2, 2);*/
ddg = BOLTZ*temp*sigma_k/k;
printf("Fitting paramaters chi^2 = %10g, Quality of fit = %10g\n",
chi2, Q);
static void normalizeACF(real *ct, real *gt, int nhb, int len)
{
- real ct_fac, gt_fac;
+ real ct_fac, gt_fac = 0;
int i;
/* Xu and Berne use the same normalization constant */
ct_fac = 1.0/ct[0];
- gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
+ if (nhb != 0)
+ {
+ gt_fac = 1.0/(real)nhb;
+ }
printf("Normalization for c(t) = %g for gh(t) = %g\n", ct_fac, gt_fac);
for (i = 0; i < len; i++)
fprintf(fp, "%10g %10g %10g %10g %10g\n",
hb->time[j]-hb->time[0], ct[j], cct[j], ght[j], kt[j]);
}
- ffclose(fp);
+ gmx_ffclose(fp);
analyse_corr(nn, hb->time, ct, ght, kt, NULL, NULL, NULL,
fit_start, temp, smooth_tail_start, oenv);
fp = opt2FILE("-hbn", nfile, fnm, "w");
if (opt2bSet("-g", nfile, fnm))
{
- fplog = ffopen(opt2fn("-g", nfile, fnm), "w");
+ fplog = gmx_ffopen(opt2fn("-g", nfile, fnm), "w");
fprintf(fplog, "# %10s %12s %12s\n", "Donor", "Hydrogen", "Acceptor");
}
else
for (i = 0; i < isize[grp]; i++)
{
fprintf(fp, (i%15) ? " " : "\n");
- fprintf(fp, " %4u", index[grp][i]+1);
+ fprintf(fp, " %4d", index[grp][i]+1);
}
fprintf(fp, "\n");
/*
{
for (j = 0; (j < hb->d.nhydro[i]); j++)
{
- fprintf(fp, " %4u %4u", hb->d.don[i]+1,
+ fprintf(fp, " %4d %4d", hb->d.don[i]+1,
hb->d.hydro[i][j]+1);
}
fprintf(fp, "\n");
if (hb->a.grp[i] == grp)
{
fprintf(fp, (i%15 && !first) ? " " : "\n");
- fprintf(fp, " %4u", hb->a.acc[i]+1);
+ fprintf(fp, " %4d", hb->a.acc[i]+1);
first = FALSE;
}
}
sprintf(as, "%s", mkatomname(atoms, aaa));
if (bContact)
{
- fprintf(fp, " %6u %6u\n", ddd+1, aaa+1);
+ fprintf(fp, " %6d %6d\n", ddd+1, aaa+1);
if (fplog)
{
fprintf(fplog, "%12s %12s\n", ds, as);
{
hhh = hb->d.hydro[i][m];
sprintf(hs, "%s", mkatomname(atoms, hhh));
- fprintf(fp, " %6u %6u %6u\n", ddd+1, hhh+1, aaa+1);
+ fprintf(fp, " %6d %6d %6d\n", ddd+1, hhh+1, aaa+1);
if (fplog)
{
fprintf(fplog, "%12s %12s %12s\n", ds, hs, as);
}
}
}
- ffclose(fp);
+ gmx_ffclose(fp);
if (fplog)
{
- ffclose(fplog);
+ gmx_ffclose(fplog);
}
}
"one particle" },
{ "-fitstart", FALSE, etREAL, {&fit_start},
"Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. With [TT]-gemfit[tt] we suggest [TT]-fitstart 0[tt]" },
- { "-fitstart", FALSE, etREAL, {&fit_start},
+ { "-fitend", FALSE, etREAL, {&fit_end},
"Time (ps) to which to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation (only with [TT]-gemfit[tt])" },
{ "-temp", FALSE, etREAL, {&temp},
"Temperature (K) for computing the Gibbs energy corresponding to HB breaking and reforming" },
{ "-merge", FALSE, etBOOL, {&bMerge},
"H-bonds between the same donor and acceptor, but with different hydrogen are treated as a single H-bond. Mainly important for the ACF." },
{ "-geminate", FALSE, etENUM, {gemType},
- "Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
+ "HIDDENUse reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
{ "-diff", FALSE, etREAL, {&D},
- "Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
+ "HIDDENDffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
#ifdef GMX_OPENMP
{ "-nthreads", FALSE, etINT, {&nThreads},
- "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
+ "Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of cores (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
#endif
};
const char *bugs[] = {
};
t_filenm fnm[] = {
{ efTRX, "-f", NULL, ffREAD },
- { efTPX, NULL, NULL, ffREAD },
+ { efTPR, NULL, NULL, ffREAD },
{ efNDX, NULL, NULL, ffOPTRD },
/* { efNDX, "-sel", "select", ffOPTRD },*/
{ efXVG, "-num", "hbnum", ffWRITE },
npargs = asize(pa);
ppa = add_acf_pargs(&npargs, pa);
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE, NFILE, fnm, npargs,
+ if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT, NFILE, fnm, npargs,
ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
{
return 0;
if (bGem)
{
printf("Geminate recombination: %s\n", gemType[gemmode]);
-#ifndef HAVE_LIBGSL
- printf("Note that some aspects of reversible geminate recombination won't work without gsl.\n");
-#endif
if (bContact)
{
if (gemmode != gemDD)
hb = mk_hbdata(bHBmap, opt2bSet("-dan", NFILE, fnm), bMerge || bContact, bGem, gemmode);
/* get topology */
- read_tpx_top(ftp2fn(efTPX, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
+ read_tpx_top(ftp2fn(efTPR, NFILE, fnm), &ir, box, &natoms, NULL, NULL, NULL, &top);
snew(grpnames, grNR);
snew(index, grNR);
close_trj(status);
if (fpnhb)
{
- ffclose(fpnhb);
+ gmx_ffclose(fpnhb);
}
/* Compute maximum possible number of different hbonds */
aver_nhb += hb->nhb[i];
aver_dist += hb->ndist[i];
}
- ffclose(fp);
+ gmx_ffclose(fp);
aver_nhb /= nframes;
aver_dist /= nframes;
/* Print HB distance distribution */
{
fprintf(fp, "%10g %10g\n", (i+0.5)*rbin, rdist[i]/(rbin*(real)sum));
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
/* Print HB angle distribution */
{
fprintf(fp, "%10g %10g\n", (i+0.5)*abin, adist[i]/(abin*(real)sum));
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
/* Print HB in alpha-helix */
}
fprintf(fp, "\n");
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
if (!bNN)
{
gmx_fatal(FARGS, "Could not initiate t_gemParams params.");
}
}
- gemstring = strdup(gemType[hb->per->gemtype]);
+ gemstring = gmx_strdup(gemType[hb->per->gemtype]);
do_hbac(opt2fn("-ac", NFILE, fnm), hb, nDump,
bMerge, bContact, fit_start, temp, r2cut > 0, smooth_tail_start, oenv,
gemstring, nThreads, NN, bBallistic, bGemFit);
{
t_matrix mat;
int id, ia, hh, x, y;
+ mat.flags = mat.y0 = 0;
if ((nframes > 0) && (hb->nrhb > 0))
{
}
fp = opt2FILE("-hbm", NFILE, fnm, "w");
write_xpm_m(fp, mat);
- ffclose(fp);
+ gmx_ffclose(fp);
for (x = 0; x < mat.nx; x++)
{
sfree(mat.matrix[x]);
if (USE_THIS_GROUP(j) )
{
sprintf(buf, "Donors %s", grpnames[j]);
- legnames[i++] = strdup(buf);
+ legnames[i++] = gmx_strdup(buf);
sprintf(buf, "Acceptors %s", grpnames[j]);
- legnames[i++] = strdup(buf);
+ legnames[i++] = gmx_strdup(buf);
}
}
if (i != nleg)
}
fprintf(fp, "\n");
}
- ffclose(fp);
+ gmx_ffclose(fp);
}
return 0;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff