From b5daddfbac614cc00cb467c515317913784d3b4b Mon Sep 17 00:00:00 2001 From: David van der Spoel Date: Fri, 7 Feb 2014 10:45:37 +0100 Subject: [PATCH] Removed possibility of compiling with GSL. Cleanup up of programs that are not functioning well anyway will be done in conjunction with a move to the new analysis framework. Removed g_kinetics. Refs #1472 Change-Id: Ic948e239fd77d8ae519a6462c1e25f908a5de014 --- CMakeLists.txt | 11 - admin/installguide/installguide.tex | 2 - cmake/FindGSL.cmake | 72 -- src/config.h.cmakein | 3 - src/gromacs/CMakeLists.txt | 2 +- src/gromacs/gmxana/geminate.c | 498 +---------- src/gromacs/gmxana/gmx_ana.h | 3 - src/gromacs/gmxana/gmx_hbond.c | 164 +--- src/gromacs/gmxana/gmx_kinetics.c | 1081 ------------------------ src/gromacs/libgromacs.pc.cmakein | 2 +- src/programs/CreateLinks.cmake.cmakein | 1 - src/programs/legacymodules.cpp | 3 - 12 files changed, 16 insertions(+), 1826 deletions(-) delete mode 100644 cmake/FindGSL.cmake delete mode 100644 src/gromacs/gmxana/gmx_kinetics.c diff --git a/CMakeLists.txt b/CMakeLists.txt index ff512eeab0..4922e03f47 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -429,17 +429,6 @@ mark_as_advanced(GMX_XML) # set(XML_LIBRARIES ${LIBXML2_LIBRARIES}) #endif() -option(GMX_GSL "Add support for gsl" OFF) -if (GMX_GSL) - find_package(GSL) - set(PKG_GSL "") - if(GSL_FOUND) - include_directories(${GSL_INCLUDE_DIR}) - set(PKG_GSL gsl) - set(HAVE_LIBGSL 1) - endif() -endif() - option(GMX_EXTRAE "Add support for tracing using EXTRAE" OFF) mark_as_advanced(GMX_EXTRAE) diff --git a/admin/installguide/installguide.tex b/admin/installguide/installguide.tex index cfee6c537c..ab7ef8753e 100644 --- a/admin/installguide/installguide.tex +++ b/admin/installguide/installguide.tex @@ -308,8 +308,6 @@ Intel's \mkl{} documentation for your system. to use third-party software for visualization, such as \vmd{} \url{http://www.ks.uiuc.edu/Research/vmd/} or \pymol{} \url{http://www.pymol.org/}. -\item A few \gromacs{} tools get some extra functionality when linked with the -GNU scientific library GSL. \item Running the \gromacs{} test suite requires \libxmltwo{} \item Building the \gromacs{} manual requires ImageMagick, pdflatex and bibtex diff --git a/cmake/FindGSL.cmake b/cmake/FindGSL.cmake deleted file mode 100644 index 26440eb647..0000000000 --- a/cmake/FindGSL.cmake +++ /dev/null @@ -1,72 +0,0 @@ -# -# This file is part of the GROMACS molecular simulation package. -# -# Copyright (c) 2009-2011, by the VOTCA Development Team (http://www.votca.org). -# Copyright (c) 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. - -# - Find gsl -# Find the native GSL headers and libraries. -# -# GSL_INCLUDE_DIRS - where to find gsl/gsl_linalg.h, etc. -# GSL_LIBRARIES - List of libraries when using gsl. -# GSL_FOUND - True if gsl found. -# - -find_package(PkgConfig QUIET) - -pkg_check_modules(PC_GSL gsl) -find_path(GSL_INCLUDE_DIR gsl/gsl_linalg.h HINTS ${PC_GSL_INCLUDE_DIRS}) -find_library(GSL_LIBRARY NAMES gsl HINTS ${PC_GSL_LIBRARY_DIRS} ) -find_library(GSLCBLAS_LIBRARY NAMES gslcblas cblas HINTS ${PC_GSL_LIBRARY_DIRS} ) - -include(FindPackageHandleStandardArgs) -# handle the QUIETLY and REQUIRED arguments and set GSL_FOUND to TRUE -# if all listed variables are TRUE -find_package_handle_standard_args(GSL DEFAULT_MSG GSL_LIBRARY GSL_INCLUDE_DIR GSLCBLAS_LIBRARY) - -set(GSL_LIBRARIES "${GSL_LIBRARY};${GSLCBLAS_LIBRARY}") -set(GSL_INCLUDE_DIRS ${GSL_INCLUDE_DIR} ) - -if (GSL_FOUND) - include(CheckLibraryExists) - check_library_exists("${GSLCBLAS_LIBRARY};${MATH_LIBRARIES}" cblas_dsyrk "" FOUND_DSYRK) - if(NOT FOUND_DSYRK) - message(FATAL_ERROR "Could not find cblas_dsyrk in ${GSLCBLAS_LIBRARY}, take a look at the error message in ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log to find out what was going wrong. If you are using a static lib (.a) make sure you have specified all dependencies of libcblas in GSLCBLAS_LIBRARY by hand (e.g. -DGSLCBLAS_LIBRARY='/path/to/libcblas.so;/path/to/libm.so')! If your gsl was build against an different version of cblas, specify it in GSLCBLAS_LIBRARY") - endif(NOT FOUND_DSYRK) - #adding MATH_LIBRARIES here to allow static libs, this does not harm us as we are anyway using it - check_library_exists("${GSL_LIBRARIES};${MATH_LIBRARIES}" gsl_linalg_QR_decomp "" FOUND_QR_DECOMP) - if(NOT FOUND_QR_DECOMP) - message(FATAL_ERROR "Could not find gsl_linalg_QR_decompx in ${GSL_LIBRARY}, take a look at the error message in ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log to find out what was going wrong. If you are using a static lib (.a) make sure you have specified all dependencies of libgsl in GSL_LIBRARY by hand (e.g. -DGSL_LIBRARY='/path/to/libgsl.so;/path/to/libm.so') !") - endif(NOT FOUND_QR_DECOMP) -endif (GSL_FOUND) - -mark_as_advanced(GSL_INCLUDE_DIR GSL_LIBRARY) diff --git a/src/config.h.cmakein b/src/config.h.cmakein index 0a75d07ac3..c0fe73d85e 100644 --- a/src/config.h.cmakein +++ b/src/config.h.cmakein @@ -223,9 +223,6 @@ /* Define to 1 if _fseeki64 (and presumably _fseeki64) exists and is declared. */ #cmakedefine HAVE__FSEEKI64 -/* Define to 1 if you have the gsl library (-lgsl). */ -#cmakedefine HAVE_LIBGSL - /* Have io.h (windows)*/ #cmakedefine HAVE_IO_H diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt index 21ba8166a6..32194e942d 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -158,7 +158,7 @@ target_link_libraries(libgromacs ${GMX_GPU_LIBRARIES} ${GMX_TNG_LIBRARIES} ${EXTRAE_LIBRARIES} ${FFT_LIBRARIES} ${LINEAR_ALGEBRA_LIBRARIES} - ${XML_LIBRARIES} ${GSL_LIBRARIES} + ${XML_LIBRARIES} ${THREAD_LIB} ${GMX_SHARED_LINKER_FLAGS}) set_target_properties(libgromacs PROPERTIES OUTPUT_NAME "gromacs${GMX_LIBS_SUFFIX}" diff --git a/src/gromacs/gmxana/geminate.c b/src/gromacs/gmxana/geminate.c index c9be94d374..bfb7153fbe 100644 --- a/src/gromacs/gmxana/geminate.c +++ b/src/gromacs/gmxana/geminate.c @@ -36,17 +36,6 @@ #include #endif -#ifdef HAVE_LIBGSL -#include -#include -#include -#include -#include -#include -#include -#include -#endif - #include #include "typedefs.h" #include "gromacs/utility/smalloc.h" @@ -55,6 +44,12 @@ #include "gromacs/utility/gmxomp.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). @@ -589,63 +584,6 @@ extern t_gemParams *init_gemParams(const double sigma, const double D, * 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. */ -#ifdef HAVE_LIBGSL -static double gemFunc_residual2(const gsl_vector *p, void *data) -{ - gemFitData *GD; - int i, iLog, nLin, nFitPoints, nData; - double r, residual2, *ctTheory, *y; - - GD = (gemFitData *)data; - nLin = GD->params->nLin; - nFitPoints = GD->params->nFitPoints; - nData = GD->nData; - residual2 = 0; - ctTheory = GD->ctTheory; - y = GD->y; - - /* Now, we'd save a lot of time by not calculating eq10v2 for all - * time points, but only those that we sample to calculate the mse. - */ - - eq10v2(GD->ctTheory, GD->doubleLogTime /* GD->time */, nFitPoints /* GD->nData */, - gsl_vector_get(p, 0), gsl_vector_get(p, 1), - GD->params); - - fixGemACF(GD->ctTheory, nFitPoints); - - /* Removing a bunch of points from the log-part. */ -#pragma omp parallel for schedule(dynamic) \ - firstprivate(nData, ctTheory, y, nFitPoints) \ - private (i, iLog, r) \ - reduction(+:residual2) \ - default(shared) - for (i = 0; i < nFitPoints; i++) - { - iLog = GD->logtime[i]; - r = log(ctTheory[i /* iLog */]); - residual2 += sqr(r-log(y[iLog])); - } - residual2 /= nFitPoints; /* Not really necessary for the fitting, but gives more meaning to the returned data. */ - /* printf("ka=%3.5f\tkd=%3.5f\trmse=%3.5f\n", gsl_vector_get(p,0), gsl_vector_get(p,1), residual2); */ - return residual2; - -/* for (i=0; iparams); */ -/* if (iLog >= nData) */ -/* gmx_fatal(FARGS, "log index out of bounds: %i", iLog); */ -/* r = ctTheory[iLog]; */ -/* residual2 += sqr(r-y[iLog]); */ - -/* } */ -/* residual2 /= GD->n; */ -/* return residual2; */ -} -#endif static real* d2r(const double *d, const int nn); @@ -662,253 +600,11 @@ extern real fitGemRecomb(double gmx_unused *ct, gemFitData *GD; char *dumpstr, dumpname[128]; - /* nmsimplex2 had convergence problems prior to gsl v1.14, - * but it's O(N) instead of O(N) operations, so let's use it if v >= 1.14 */ -#ifdef HAVE_LIBGSL - gsl_multimin_fminimizer *s; - gsl_vector *x, *dx; /* parameters and initial step size */ - gsl_multimin_function fitFunc; -#ifdef GSL_MAJOR_VERSION -#ifdef GSL_MINOR_VERSION -#if ((GSL_MAJOR_VERSION == 1 && GSL_MINOR_VERSION >= 14) || \ - (GSL_MAJOR_VERSION > 1)) - const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2; -#else - const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; -#endif /* #if ... */ -#endif /* GSL_MINOR_VERSION */ -#else - const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; -#endif /* GSL_MAJOR_VERSION */ - fprintf(stdout, "Will fit ka and kd to the ACF according to the reversible geminate recombination model.\n"); -#else /* HAVE_LIBGSL */ - fprintf(stderr, "Sorry, can't do reversible geminate recombination without gsl. " - "Recompile using --with-gsl.\n"); + missing_code_message(); return -1; -#endif /* HAVE_LIBGSL */ - -#ifdef HAVE_LIBGSL -#ifdef GMX_OPENMP - nThreads = gmx_omp_get_max_threads(); - gmx_omp_set_num_threads(nThreads); - fprintf(stdout, "We will be using %i threads.\n", nThreads); -#endif - - iter = 0; - status = 0; - maxiter = 100; - tol = 1e-10; - - p = 2; /* Number of parameters to fit. ka and kd. */ - n = params->nFitPoints; /* params->nLin*2 */; /* Number of points in the reduced dataset */ - - if (params->D <= 0) - { - fprintf(stderr, "Fitting of D is not implemented yet. It must be provided on the command line.\n"); - return -1; - } - -/* if (nDatan = n; - GD->y = ct; - GD->ctTheory = NULL; - snew(GD->ctTheory, nData); - GD->LinLog = NULL; - snew(GD->LinLog, n); - GD->time = time; - GD->ka = 0; - GD->kd = 0; - GD->tDelta = time[1]-time[0]; - GD->nData = nData; - GD->params = params; - snew(GD->logtime, params->nFitPoints); - snew(GD->doubleLogTime, params->nFitPoints); - - for (i = 0; i < params->nFitPoints; i++) - { - GD->doubleLogTime[i] = (double)(getLogIndex(i, params)); - GD->logtime[i] = (int)(GD->doubleLogTime[i]); - GD->doubleLogTime[i] *= GD->tDelta; - - if (GD->logtime[i] >= nData) - { - fprintf(stderr, "Ayay. It seems we're indexing out of bounds.\n"); - params->nFitPoints = i; - } - } - - fitFunc.f = &gemFunc_residual2; - fitFunc.n = 2; - fitFunc.params = (void*)GD; - - x = gsl_vector_alloc (fitFunc.n); - dx = gsl_vector_alloc (fitFunc.n); - gsl_vector_set (x, 0, 25); - gsl_vector_set (x, 1, 0.5); - gsl_vector_set (dx, 0, 0.1); - gsl_vector_set (dx, 1, 0.01); - - - s = gsl_multimin_fminimizer_alloc (T, fitFunc.n); - gsl_multimin_fminimizer_set (s, &fitFunc, x, dx); - gsl_vector_free (x); - gsl_vector_free (dx); - - do - { - iter++; - status = gsl_multimin_fminimizer_iterate (s); - - if (status != 0) - { - gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s:\n \"%s\"\n", - gsl_multimin_fminimizer_name(s), gsl_strerror(status)); - } - - d2 = gsl_multimin_fminimizer_minimum(s); - size = gsl_multimin_fminimizer_size(s); - params->ka = gsl_vector_get (s->x, 0); - params->kd = gsl_vector_get (s->x, 1); - - if (status) - { - fprintf(stderr, "%s\n", gsl_strerror(status)); - break; - } - - status = gsl_multimin_test_size(size, tol); - - if (status == GSL_SUCCESS) - { - fprintf(stdout, "Converged to minimum at\n"); - } - - printf ("iter %5d: ka = %2.5f kd = %2.5f f() = %7.3f size = %.3f chi2 = %2.5f\n", - iter, - params->ka, - params->kd, - s->fval, size, d2); - - if (iter%1 == 0) - { - eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params); - /* fixGemACF(GD->ctTheory, nFitPoints); */ - sprintf(dumpname, "Iter_%i.xvg", iter); - for (i = 0; i < GD->nData; i++) - { - dumpdata[i] = (real)(GD->ctTheory[i]); - if (!gmx_isfinite(dumpdata[i])) - { - gmx_fatal(FARGS, "Non-finite value in acf."); - } - } - dumpN(dumpdata, GD->nData, dumpname); - } - } - while ((status == GSL_CONTINUE) && (iter < maxiter)); - - /* /\* Calculate the theoretical ACF from the parameters one last time. *\/ */ - eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params); - *ctFit = GD->ctTheory; - sfree(GD); - gsl_multimin_fminimizer_free (s); - - - return d2; - -#endif /* HAVE_LIBGSL */ -} - -#ifdef HAVE_LIBGSL -static int balFunc_f(const gsl_vector *x, void *data, gsl_vector *f) -{ - /* C + sum{ A_i * exp(-B_i * t) }*/ - - balData *BD; - int n, i, j, nexp; - double *y, *A, *B, C, /* There are the parameters to be optimized. */ - t, ct; - - BD = (balData *)data; - n = BD->n; - nexp = BD->nexp; - y = BD->y, - snew(A, nexp); - snew(B, nexp); - - for (i = 0; i < nexp; i++) - { - A[i] = gsl_vector_get(x, i*2); - B[i] = gsl_vector_get(x, i*2+1); - } - C = gsl_vector_get(x, nexp*2); - - for (i = 0; i < n; i++) - { - t = i*BD->tDelta; - ct = 0; - for (j = 0; j < nexp; j++) - { - ct += A[j] * exp(B[j] * t); - } - ct += C; - gsl_vector_set (f, i, ct - y[i]); - } - return GSL_SUCCESS; } -/* The derivative stored in jacobian form (J)*/ -static int balFunc_df(const gsl_vector *params, void *data, gsl_matrix *J) -{ - balData *BD; - size_t n, i, j; - double *y, *A, *B, C, /* There are the parameters. */ - t; - int nexp; - - BD = (balData*)data; - n = BD->n; - y = BD->y; - nexp = BD->nexp; - - snew(A, nexp); - snew(B, nexp); - - for (i = 0; i < nexp; i++) - { - A[i] = gsl_vector_get(params, i*2); - B[i] = gsl_vector_get(params, i*2+1); - } - C = gsl_vector_get(params, nexp*2); - for (i = 0; i < n; i++) - { - t = i*BD->tDelta; - for (j = 0; j < nexp; j++) - { - gsl_matrix_set (J, i, j*2, exp(B[j]*t)); /* df(t)/dA_j */ - gsl_matrix_set (J, i, j*2+1, A[j]*t*exp(B[j]*t)); /* df(t)/dB_j */ - } - gsl_matrix_set (J, i, nexp*2, 1); /* df(t)/dC */ - } - return GSL_SUCCESS; -} - -/* Calculation of the function and its derivative */ -static int balFunc_fdf(const gsl_vector *params, void *data, - gsl_vector *f, gsl_matrix *J) -{ - balFunc_f(params, data, f); - balFunc_df(params, data, J); - return GSL_SUCCESS; -} -#endif /* HAVE_LIBGSL */ /* Removes the ballistic term from the beginning of the ACF, * just like in Omer's paper. @@ -916,8 +612,7 @@ static int balFunc_fdf(const gsl_vector *params, void *data, extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tMax, int nexp, gmx_bool gmx_unused bDerivative) { - /* Use nonlinear regression with GSL instead. - * Fit with 4 exponentials and one constant term, + /* Fit with 4 exponentials and one constant term, * subtract the fatest exponential. */ int nData, i, status, iter; @@ -939,183 +634,8 @@ extern void takeAwayBallistic(double gmx_unused *ct, double *t, int len, real tM p = nexp*2+1; /* Number of parameters. */ -#ifdef HAVE_LIBGSL - const gsl_multifit_fdfsolver_type *T - = gsl_multifit_fdfsolver_lmsder; - - gsl_multifit_fdfsolver *s; /* The solver itself. */ - gsl_multifit_function_fdf fitFunction; /* The function to be fitted. */ - gsl_matrix *covar; /* Covariance matrix for the parameters. - * We'll not use the result, though. */ - gsl_vector_view theParams; - - nData = 0; - do - { - nData++; - } - while (t[nData] < tMax+t[0] && nData < len); - - guess = NULL; - n = nData; - - snew(guess, p); - snew(A, p); - covar = gsl_matrix_alloc (p, p); - - /* Set up an initial gess for the parameters. - * The solver is somewhat sensitive to the initial guess, - * but this worked fine for a TIP3P box with -geminate dd - * EDIT: In fact, this seems like a good starting pont for other watermodels too. */ - for (i = 0; i < nexp; i++) - { - guess[i*2] = 0.1; - guess[i*2+1] = -0.5 + (((double)i)/nexp - 0.5)*0.3; - } - guess[nexp * 2] = 0.01; - - theParams = gsl_vector_view_array(guess, p); - - snew(BD, 1); - BD->n = n; - BD->y = ct; - BD->tDelta = t[1]-t[0]; - BD->nexp = nexp; - - fitFunction.f = &balFunc_f; - fitFunction.df = &balFunc_df; - fitFunction.fdf = &balFunc_fdf; - fitFunction.n = nData; - fitFunction.p = p; - fitFunction.params = BD; - - s = gsl_multifit_fdfsolver_alloc (T, nData, p); - if (s == NULL) - { - gmx_fatal(FARGS, "Could not set up the nonlinear solver."); - } - - gsl_multifit_fdfsolver_set(s, &fitFunction, &theParams.vector); - - /* \=============================================/ */ - - iter = 0; - do - { - iter++; - status = gsl_multifit_fdfsolver_iterate (s); - - if (status) - { - break; - } - status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4); - } - while (iter < 5000 && status == GSL_CONTINUE); - - if (iter == 5000) - { - fprintf(stderr, "The non-linear fitting did not converge in 5000 steps.\n" - "Check the quality of the fit!\n"); - } - else - { - fprintf(stderr, "Non-linear fitting of ballistic term converged in %d steps.\n\n", (int)iter); - } - for (i = 0; i < nexp; i++) - { - fprintf(stdout, "%c * exp(%c * t) + ", 'A'+(char)i*2, 'B'+(char)i*2); - } - - fprintf(stdout, "%c\n", 'A'+(char)nexp*2); - fprintf(stdout, "Here are the actual numbers for A-%c:\n", 'A'+nexp*2); - - for (i = 0; i < nexp; i++) - { - A[i*2] = gsl_vector_get(s->x, i*2); - A[i*2+1] = gsl_vector_get(s->x, i*2+1); - fprintf(stdout, " %g*exp(%g * x) +", A[i*2], A[i*2+1]); - } - A[i*2] = gsl_vector_get(s->x, i*2); /* The last and constant term */ - fprintf(stdout, " %g\n", A[i*2]); - - fflush(stdout); - - /* Implement some check for parameter quality */ - for (i = 0; i < nexp; i++) - { - if (A[i*2] < 0 || A[i*2] > 1) - { - fprintf(stderr, "WARNING: ----------------------------------\n" - " | A coefficient does not lie within [0,1].\n" - " | This may or may not be a problem.\n" - " | Double check the quality of the fit!\n"); - } - if (A[i*2+1] > 0) - { - fprintf(stderr, "WARNING: ----------------------------------\n" - " | One factor in the exponent is positive.\n" - " | This could be a problem if the coefficient\n" - " | is large. Double check the quality of the fit!\n"); - } - } - if (A[i*2] < 0 || A[i*2] > 1) - { - fprintf(stderr, "WARNING: ----------------------------------\n" - " | The constant term does not lie within [0,1].\n" - " | This may or may not be a problem.\n" - " | Double check the quality of the fit!\n"); - } - - /* Sort the terms */ - sorted = (nexp > 1) ? FALSE : TRUE; - while (!sorted) - { - sorted = TRUE; - for (i = 0; i < nexp-1; i++) - { - ddt[0] = A[i*2] * A[i*2+1]; - ddt[1] = A[i*2+2] * A[i*2+3]; - - if ((bDerivative && (ddt[0] < 0 && ddt[1] < 0 && ddt[0] > ddt[1])) || /* Compare derivative at t=0... */ - (!bDerivative && (A[i*2+1] > A[i*2+3]))) /* Or just the coefficient in the exponent */ - { - sorted = FALSE; - a[0] = A[i*2]; /* coefficient */ - a[1] = A[i*2+1]; /* parameter in the exponent */ - - A[i*2] = A[i*2+2]; - A[i*2+1] = A[i*2+3]; - - A[i*2+2] = a[0]; - A[i*2+3] = a[1]; - } - } - } - - /* Subtract the fastest component */ - fprintf(stdout, "Fastest component is %g * exp(%g * t)\n" - "Subtracting fastest component from ACF.\n", A[0], A[1]); - - for (i = 0; i < len; i++) - { - ct[i] = (ct[i] - A[0] * exp(A[1] * i*BD->tDelta)) / (1-A[0]); - } - - sfree(guess); - sfree(A); - - gsl_multifit_fdfsolver_free(s); - gsl_matrix_free(covar); - fflush(stdout); - -#else - /* We have no gsl. */ - fprintf(stderr, "Sorry, can't take away ballistic component without gsl. " - "Recompile using --with-gsl.\n"); + missing_code_message(); return; -#endif /* HAVE_LIBGSL */ - } extern void dumpN(const real *e, const int nn, char *fn) diff --git a/src/gromacs/gmxana/gmx_ana.h b/src/gromacs/gmxana/gmx_ana.h index 1b2e58bc26..acac0e9beb 100644 --- a/src/gromacs/gmxana/gmx_ana.h +++ b/src/gromacs/gmxana/gmx_ana.h @@ -147,9 +147,6 @@ gmx_helixorient(int argc, char *argv[]); int gmx_hydorder(int argc, char *argv[]); -int -gmx_kinetics(int argc, char *argv[]); - int gmx_make_edi(int argc, char *argv[]); diff --git a/src/gromacs/gmxana/gmx_hbond.c b/src/gromacs/gmxana/gmx_hbond.c index 5c0206eeb6..a9b8942786 100644 --- a/src/gromacs/gmxana/gmx_hbond.c +++ b/src/gromacs/gmxana/gmx_hbond.c @@ -2317,157 +2317,6 @@ typedef struct { real *t, *ct, *nt, *kt, *sigma_ct, *sigma_nt, *sigma_kt; } t_luzar; -#ifdef HAVE_LIBGSL -#include -#include -#include - -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, @@ -2503,13 +2352,13 @@ static real compute_weighted_rates(int n, real t[], real ct[], real nt[], 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]); @@ -2597,7 +2446,7 @@ void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], 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); @@ -3772,9 +3621,9 @@ int gmx_hbond(int argc, char *argv[]) { "-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)"}, @@ -3902,9 +3751,6 @@ int gmx_hbond(int argc, char *argv[]) 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) diff --git a/src/gromacs/gmxana/gmx_kinetics.c b/src/gromacs/gmxana/gmx_kinetics.c deleted file mode 100644 index 2854ded641..0000000000 --- a/src/gromacs/gmxana/gmx_kinetics.c +++ /dev/null @@ -1,1081 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * 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. - * - * 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. - */ -#ifdef HAVE_CONFIG_H -#include -#endif -#include -#include -#include -#include "gromacs/commandline/pargs.h" -#include "sysstuff.h" -#include "typedefs.h" -#include "gromacs/utility/smalloc.h" -#include "macros.h" -#include "gmx_fatal.h" -#include "vec.h" -#include "copyrite.h" -#include "gromacs/fileio/futil.h" -#include "readinp.h" -#include "txtdump.h" -#include "gstat.h" -#include "xvgr.h" -#include "physics.h" -#include "gmx_ana.h" - -enum { - epAuf, epEuf, epAfu, epEfu, epNR -}; -enum { - eqAif, eqEif, eqAfi, eqEfi, eqAui, eqEui, eqAiu, eqEiu, eqNR -}; -static char *eep[epNR] = { "Af", "Ef", "Au", "Eu" }; -static char *eeq[eqNR] = { "Aif", "Eif", "Afi", "Efi", "Aui", "Eui", "Aiu", "Eiu" }; - -typedef struct { - int nreplica; /* Number of replicas in the calculation */ - int nframe; /* Number of time frames */ - int nstate; /* Number of states the system can be in, e.g. F,I,U */ - int nparams; /* Is 2, 4 or 8 */ - gmx_bool *bMask; /* Determine whether this replica is part of the d2 comp. */ - gmx_bool bSum; - gmx_bool bDiscrete; /* Use either discrete folding (0/1) or a continuous */ - /* criterion */ - int nmask; /* Number of replicas taken into account */ - real dt; /* Timestep between frames */ - int j0, j1; /* Range of frames used in calculating delta */ - real **temp, **data, **data2; - int **state; /* State index running from 0 (F) to nstate-1 (U) */ - real **beta, **fcalt, **icalt; - real *time, *sumft, *sumit, *sumfct, *sumict; - real *params; - real *d2_replica; -} t_remd_data; - -#ifdef HAVE_LIBGSL -#include - -static char *itoa(int i) -{ - static char ptr[12]; - - sprintf(ptr, "%d", i); - return ptr; -} - -static char *epnm(int nparams, int index) -{ - static char buf[32], from[8], to[8]; - int nn, ni, ii; - - range_check(index, 0, nparams); - if ((nparams == 2) || (nparams == 4)) - { - return eep[index]; - } - else if ((nparams > 4) && (nparams % 4 == 0)) - { - return eeq[index]; - } - else - { - gmx_fatal(FARGS, "Don't know how to handle %d parameters", nparams); - } - - return NULL; -} - -static gmx_bool bBack(t_remd_data *d) -{ - return (d->nparams > 2); -} - -static real is_folded(t_remd_data *d, int irep, int jframe) -{ - if (d->state[irep][jframe] == 0) - { - return 1.0; - } - else - { - return 0.0; - } -} - -static real is_unfolded(t_remd_data *d, int irep, int jframe) -{ - if (d->state[irep][jframe] == d->nstate-1) - { - return 1.0; - } - else - { - return 0.0; - } -} - -static real is_intermediate(t_remd_data *d, int irep, int jframe) -{ - if ((d->state[irep][jframe] == 1) && (d->nstate > 2)) - { - return 1.0; - } - else - { - return 0.0; - } -} - -static void integrate_dfdt(t_remd_data *d) -{ - int i, j; - double beta, ddf, ddi, df, db, fac, sumf, sumi, area; - - d->sumfct[0] = 0; - d->sumict[0] = 0; - for (i = 0; (i < d->nreplica); i++) - { - if (d->bMask[i]) - { - if (d->bDiscrete) - { - ddf = 0.5*d->dt*is_folded(d, i, 0); - ddi = 0.5*d->dt*is_intermediate(d, i, 0); - } - else - { - ddf = 0.5*d->dt*d->state[i][0]; - ddi = 0.0; - } - d->fcalt[i][0] = ddf; - d->icalt[i][0] = ddi; - d->sumfct[0] += ddf; - d->sumict[0] += ddi; - } - } - for (j = 1; (j < d->nframe); j++) - { - if (j == d->nframe-1) - { - fac = 0.5*d->dt; - } - else - { - fac = d->dt; - } - sumf = sumi = 0; - for (i = 0; (i < d->nreplica); i++) - { - if (d->bMask[i]) - { - beta = d->beta[i][j]; - if ((d->nstate <= 2) || d->bDiscrete) - { - if (d->bDiscrete) - { - df = (d->params[epAuf]*exp(-beta*d->params[epEuf])* - is_unfolded(d, i, j)); - } - else - { - area = (d->data2 ? d->data2[i][j] : 1.0); - df = area*d->params[epAuf]*exp(-beta*d->params[epEuf]); - } - if (bBack(d)) - { - db = 0; - if (d->bDiscrete) - { - db = (d->params[epAfu]*exp(-beta*d->params[epEfu])* - is_folded(d, i, j)); - } - else - { - gmx_fatal(FARGS, "Back reaction not implemented with continuous"); - } - ddf = fac*(df-db); - } - else - { - ddf = fac*df; - } - d->fcalt[i][j] = d->fcalt[i][j-1] + ddf; - sumf += ddf; - } - else - { - ddf = fac*((d->params[eqAif]*exp(-beta*d->params[eqEif])* - is_intermediate(d, i, j)) - - (d->params[eqAfi]*exp(-beta*d->params[eqEfi])* - is_folded(d, i, j))); - ddi = fac*((d->params[eqAui]*exp(-beta*d->params[eqEui])* - is_unfolded(d, i, j)) - - (d->params[eqAiu]*exp(-beta*d->params[eqEiu])* - is_intermediate(d, i, j))); - d->fcalt[i][j] = d->fcalt[i][j-1] + ddf; - d->icalt[i][j] = d->icalt[i][j-1] + ddi; - sumf += ddf; - sumi += ddi; - } - } - } - d->sumfct[j] = d->sumfct[j-1] + sumf; - d->sumict[j] = d->sumict[j-1] + sumi; - } - if (debug) - { - fprintf(debug, "@type xy\n"); - for (j = 0; (j < d->nframe); j++) - { - fprintf(debug, "%8.3f %12.5e\n", d->time[j], d->sumfct[j]); - } - fprintf(debug, "&\n"); - } -} - -static void sum_ft(t_remd_data *d) -{ - int i, j; - double fac; - - for (j = 0; (j < d->nframe); j++) - { - d->sumft[j] = 0; - d->sumit[j] = 0; - if ((j == 0) || (j == d->nframe-1)) - { - fac = d->dt*0.5; - } - else - { - fac = d->dt; - } - for (i = 0; (i < d->nreplica); i++) - { - if (d->bMask[i]) - { - if (d->bDiscrete) - { - d->sumft[j] += fac*is_folded(d, i, j); - d->sumit[j] += fac*is_intermediate(d, i, j); - } - else - { - d->sumft[j] += fac*d->state[i][j]; - } - } - } - } -} - -static double calc_d2(t_remd_data *d) -{ - int i, j; - double dd2, d2 = 0, dr2, tmp; - - integrate_dfdt(d); - - if (d->bSum) - { - for (j = d->j0; (j < d->j1); j++) - { - if (d->bDiscrete) - { - d2 += sqr(d->sumft[j]-d->sumfct[j]); - if (d->nstate > 2) - { - d2 += sqr(d->sumit[j]-d->sumict[j]); - } - } - else - { - d2 += sqr(d->sumft[j]-d->sumfct[j]); - } - } - } - else - { - for (i = 0; (i < d->nreplica); i++) - { - dr2 = 0; - if (d->bMask[i]) - { - for (j = d->j0; (j < d->j1); j++) - { - tmp = sqr(is_folded(d, i, j)-d->fcalt[i][j]); - d2 += tmp; - dr2 += tmp; - if (d->nstate > 2) - { - tmp = sqr(is_intermediate(d, i, j)-d->icalt[i][j]); - d2 += tmp; - dr2 += tmp; - } - } - d->d2_replica[i] = dr2/(d->j1-d->j0); - } - } - } - dd2 = (d2/(d->j1-d->j0))/(d->bDiscrete ? d->nmask : 1); - - return dd2; -} - -static double my_f(const gsl_vector *v, void *params) -{ - t_remd_data *d = (t_remd_data *) params; - double penalty = 0; - int i; - - for (i = 0; (i < d->nparams); i++) - { - d->params[i] = gsl_vector_get(v, i); - if (d->params[i] < 0) - { - penalty += 10; - } - } - if (penalty > 0) - { - return penalty; - } - else - { - return calc_d2(d); - } -} - -static void optimize_remd_parameters(t_remd_data *d, 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 = d->nparams; - my_func.params = (void *) d; - - /* Starting point */ - x = gsl_vector_alloc (my_func.n); - for (i = 0; (i < my_func.n); i++) - { - gsl_vector_set (x, i, d->params[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.1*d->params[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); - - printf ("%5s", "Iter"); - for (i = 0; (i < my_func.n); i++) - { - printf(" %12s", epnm(my_func.n, i)); - } - printf (" %12s %12s\n", "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) - { - printf ("Minimum found using %s at:\n", - gsl_multimin_fminimizer_name(s)); - } - - printf ("%5d", iter); - for (i = 0; (i < my_func.n); i++) - { - printf(" %12.4e", gsl_vector_get (s->x, i)); - } - printf (" %12.4e %12.4e\n", size, d2); - } - while ((status == GSL_CONTINUE) && (iter < maxiter)); - - gsl_multimin_fminimizer_free (s); -} - -static void preprocess_remd(FILE *fp, t_remd_data *d, real cutoff, real tref, - real ucut, gmx_bool bBack, real Euf, real Efu, - real Ei, real t0, real t1, gmx_bool bSum, gmx_bool bDiscrete, - int nmult) -{ - int i, j, ninter; - real dd, tau_f, tau_u; - - ninter = (ucut > cutoff) ? 1 : 0; - if (ninter && (ucut <= cutoff)) - { - gmx_fatal(FARGS, "You have requested an intermediate but the cutoff for intermediates %f is smaller than the normal cutoff(%f)", ucut, cutoff); - } - - if (!bBack) - { - d->nparams = 2; - d->nstate = 2; - } - else - { - d->nparams = 4*(1+ninter); - d->nstate = 2+ninter; - } - d->bSum = bSum; - d->bDiscrete = bDiscrete; - snew(d->beta, d->nreplica); - snew(d->state, d->nreplica); - snew(d->bMask, d->nreplica); - snew(d->d2_replica, d->nreplica); - snew(d->sumft, d->nframe); - snew(d->sumit, d->nframe); - snew(d->sumfct, d->nframe); - snew(d->sumict, d->nframe); - snew(d->params, d->nparams); - snew(d->fcalt, d->nreplica); - snew(d->icalt, d->nreplica); - - /* convert_times(d->nframe,d->time); */ - - if (t0 < 0) - { - d->j0 = 0; - } - else - { - for (d->j0 = 0; (d->j0 < d->nframe) && (d->time[d->j0] < t0); d->j0++) - { - ; - } - } - if (t1 < 0) - { - d->j1 = d->nframe; - } - else - { - for (d->j1 = 0; (d->j1 < d->nframe) && (d->time[d->j1] < t1); d->j1++) - { - ; - } - } - if ((d->j1-d->j0) < d->nparams+2) - { - gmx_fatal(FARGS, "Start (%f) or end time (%f) for fitting inconsistent. Reduce t0, increase t1 or supply more data", t0, t1); - } - fprintf(fp, "Will optimize from %g to %g\n", - d->time[d->j0], d->time[d->j1-1]); - d->nmask = d->nreplica; - for (i = 0; (i < d->nreplica); i++) - { - snew(d->beta[i], d->nframe); - snew(d->state[i], d->nframe); - snew(d->fcalt[i], d->nframe); - snew(d->icalt[i], d->nframe); - d->bMask[i] = TRUE; - for (j = 0; (j < d->nframe); j++) - { - d->beta[i][j] = 1.0/(BOLTZ*d->temp[i][j]); - dd = d->data[i][j]; - if (bDiscrete) - { - if (dd <= cutoff) - { - d->state[i][j] = 0; - } - else if ((ucut > cutoff) && (dd <= ucut)) - { - d->state[i][j] = 1; - } - else - { - d->state[i][j] = d->nstate-1; - } - } - else - { - d->state[i][j] = dd*nmult; - } - } - } - sum_ft(d); - - /* Assume forward rate constant is half the total time in this - * simulation and backward is ten times as long */ - if (bDiscrete) - { - tau_f = d->time[d->nframe-1]; - tau_u = 4*tau_f; - d->params[epEuf] = Euf; - d->params[epAuf] = exp(d->params[epEuf]/(BOLTZ*tref))/tau_f; - if (bBack) - { - d->params[epEfu] = Efu; - d->params[epAfu] = exp(d->params[epEfu]/(BOLTZ*tref))/tau_u; - if (ninter > 0) - { - d->params[eqEui] = Ei; - d->params[eqAui] = exp(d->params[eqEui]/(BOLTZ*tref))/tau_u; - d->params[eqEiu] = Ei; - d->params[eqAiu] = exp(d->params[eqEiu]/(BOLTZ*tref))/tau_u; - } - } - else - { - d->params[epAfu] = 0; - d->params[epEfu] = 0; - } - } - else - { - d->params[epEuf] = Euf; - if (d->data2) - { - d->params[epAuf] = 0.1; - } - else - { - d->params[epAuf] = 20.0; - } - } -} - -static real tau(real A, real E, real T) -{ - return exp(E/(BOLTZ*T))/A; -} - -static real folded_fraction(t_remd_data *d, real tref) -{ - real tauf, taub; - - tauf = tau(d->params[epAuf], d->params[epEuf], tref); - taub = tau(d->params[epAfu], d->params[epEfu], tref); - - return (taub/(tauf+taub)); -} - -static void print_tau(FILE *gp, t_remd_data *d, real tref) -{ - real tauf, taub, ddd, fff, DG, DH, TDS, Tm, Tb, Te, Fb, Fe, Fm; - int i, np = d->nparams; - - ddd = calc_d2(d); - fprintf(gp, "Final value for Chi2 = %12.5e (%d replicas)\n", ddd, d->nmask); - tauf = tau(d->params[epAuf], d->params[epEuf], tref); - fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n", - epnm(np, epAuf), d->params[epAuf], - epnm(np, epEuf), d->params[epEuf]); - if (bBack(d)) - { - taub = tau(d->params[epAfu], d->params[epEfu], tref); - fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n", - epnm(np, epAfu), d->params[epAfu], - epnm(np, epEfu), d->params[epEfu]); - fprintf(gp, "Equilibrium properties at T = %g\n", tref); - fprintf(gp, "tau_f = %8.3f ns, tau_b = %8.3f ns\n", tauf/1000, taub/1000); - fff = taub/(tauf+taub); - DG = BOLTZ*tref*log(fff/(1-fff)); - DH = d->params[epEfu]-d->params[epEuf]; - TDS = DH-DG; - fprintf(gp, "Folded fraction F = %8.3f\n", fff); - fprintf(gp, "Unfolding energies: DG = %8.3f DH = %8.3f TDS = %8.3f\n", - DG, DH, TDS); - Tb = 260; - Te = 420; - Tm = 0; - Fm = 0; - Fb = folded_fraction(d, Tb); - Fe = folded_fraction(d, Te); - while ((Te-Tb > 0.001) && (Fm != 0.5)) - { - Tm = 0.5*(Tb+Te); - Fm = folded_fraction(d, Tm); - if (Fm > 0.5) - { - Fb = Fm; - Tb = Tm; - } - else if (Fm < 0.5) - { - Te = Tm; - Fe = Fm; - } - } - if ((Fb-0.5)*(Fe-0.5) <= 0) - { - fprintf(gp, "Melting temperature Tm = %8.3f K\n", Tm); - } - else - { - fprintf(gp, "No melting temperature detected between 260 and 420K\n"); - } - if (np > 4) - { - char *ptr; - fprintf(gp, "Data for intermediates at T = %g\n", tref); - fprintf(gp, "%8s %10s %10s %10s\n", "Name", "A", "E", "tau"); - for (i = 0; (i < np/2); i++) - { - tauf = tau(d->params[2*i], d->params[2*i+1], tref); - ptr = epnm(d->nparams, 2*i); - fprintf(gp, "%8s %10.3e %10.3e %10.3e\n", ptr+1, - d->params[2*i], d->params[2*i+1], tauf/1000); - } - } - } - else - { - fprintf(gp, "Equilibrium properties at T = %g\n", tref); - fprintf(gp, "tau_f = %8.3f\n", tauf); - } -} - -static void dump_remd_parameters(FILE *gp, t_remd_data *d, const char *fn, - const char *fn2, const char *rfn, - const char *efn, const char *mfn, int skip, real tref, - output_env_t oenv) -{ - FILE *fp, *hp; - int i, j, np = d->nparams; - real rhs, tauf, taub, fff, DG; - real *params; - const char *leg[] = { "Measured", "Fit", "Difference" }; - const char *mleg[] = { "Folded fraction", "DG (kJ/mole)"}; - char **rleg; - real fac[] = { 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03 }; -#define NFAC asize(fac) - real d2[NFAC]; - double norm; - - integrate_dfdt(d); - print_tau(gp, d, tref); - norm = (d->bDiscrete ? 1.0/d->nmask : 1.0); - - if (fn) - { - fp = xvgropen(fn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv); - xvgr_legend(fp, asize(leg), leg, oenv); - for (i = 0; (i < d->nframe); i++) - { - if ((skip <= 0) || ((i % skip) == 0)) - { - fprintf(fp, "%12.5e %12.5e %12.5e %12.5e\n", d->time[i], - d->sumft[i]*norm, d->sumfct[i]*norm, - (d->sumft[i]-d->sumfct[i])*norm); - } - } - gmx_ffclose(fp); - } - if (!d->bSum && rfn) - { - snew(rleg, d->nreplica*2); - for (i = 0; (i < d->nreplica); i++) - { - snew(rleg[2*i], 32); - snew(rleg[2*i+1], 32); - sprintf(rleg[2*i], "\\f{4}F(t) %d", i); - sprintf(rleg[2*i+1], "\\f{12}F \\f{4}(t) %d", i); - } - fp = xvgropen(rfn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv); - xvgr_legend(fp, d->nreplica*2, (const char**)rleg, oenv); - for (j = 0; (j < d->nframe); j++) - { - if ((skip <= 0) || ((j % skip) == 0)) - { - fprintf(fp, "%12.5e", d->time[j]); - for (i = 0; (i < d->nreplica); i++) - { - fprintf(fp, " %5f %9.2e", is_folded(d, i, j), d->fcalt[i][j]); - } - fprintf(fp, "\n"); - } - } - gmx_ffclose(fp); - } - - if (fn2 && (d->nstate > 2)) - { - fp = xvgropen(fn2, "Optimized fit to data", "Time (ps)", - "Fraction Intermediate", oenv); - xvgr_legend(fp, asize(leg), leg, oenv); - for (i = 0; (i < d->nframe); i++) - { - if ((skip <= 0) || ((i % skip) == 0)) - { - fprintf(fp, "%12.5e %12.5e %12.5e %12.5e\n", d->time[i], - d->sumit[i]*norm, d->sumict[i]*norm, - (d->sumit[i]-d->sumict[i])*norm); - } - } - gmx_ffclose(fp); - } - if (mfn) - { - if (bBack(d)) - { - fp = xvgropen(mfn, "Melting curve", "T (K)", "", oenv); - xvgr_legend(fp, asize(mleg), mleg, oenv); - for (i = 260; (i <= 420); i++) - { - tauf = tau(d->params[epAuf], d->params[epEuf], 1.0*i); - taub = tau(d->params[epAfu], d->params[epEfu], 1.0*i); - fff = taub/(tauf+taub); - DG = BOLTZ*i*log(fff/(1-fff)); - fprintf(fp, "%5d %8.3f %8.3f\n", i, fff, DG); - } - gmx_ffclose(fp); - } - } - - if (efn) - { - snew(params, d->nparams); - for (i = 0; (i < d->nparams); i++) - { - params[i] = d->params[i]; - } - - hp = xvgropen(efn, "Chi2 as a function of relative parameter", - "Fraction", "Chi2", oenv); - for (j = 0; (j < d->nparams); j++) - { - /* Reset all parameters to optimized values */ - fprintf(hp, "@type xy\n"); - for (i = 0; (i < d->nparams); i++) - { - d->params[i] = params[i]; - } - /* Now modify one of them */ - for (i = 0; (i < NFAC); i++) - { - d->params[j] = fac[i]*params[j]; - d2[i] = calc_d2(d); - fprintf(gp, "%s = %12g d2 = %12g\n", epnm(np, j), d->params[j], d2[i]); - fprintf(hp, "%12g %12g\n", fac[i], d2[i]); - } - fprintf(hp, "&\n"); - } - gmx_ffclose(hp); - for (i = 0; (i < d->nparams); i++) - { - d->params[i] = params[i]; - } - sfree(params); - } - if (!d->bSum) - { - for (i = 0; (i < d->nreplica); i++) - { - fprintf(gp, "Chi2[%3d] = %8.2e\n", i, d->d2_replica[i]); - } - } -} -#endif /*HAVE_LIBGSL*/ - -int gmx_kinetics(int argc, char *argv[]) -{ - const char *desc[] = { - "[THISMODULE] reads two [TT].xvg[tt] files, each one containing data for N replicas.", - "The first file contains the temperature of each replica at each timestep,", - "and the second contains real values that can be interpreted as", - "an indicator for folding. If the value in the file is larger than", - "the cutoff it is taken to be unfolded and the other way around.[PAR]", - "From these data an estimate of the forward and backward rate constants", - "for folding is made at a reference temperature. In addition,", - "a theoretical melting curve and free energy as a function of temperature", - "are printed in an [TT].xvg[tt] file.[PAR]", - "The user can give a max value to be regarded as intermediate", - "([TT]-ucut[tt]), which, when given will trigger the use of an intermediate state", - "in the algorithm to be defined as those structures that have", - "cutoff < DATA < ucut. Structures with DATA values larger than ucut will", - "not be regarded as potential folders. In this case 8 parameters are optimized.[PAR]", - "The average fraction foled is printed in an [TT].xvg[tt] file together with the fit to it.", - "If an intermediate is used a further file will show the build of the intermediate and the fit to that process.[PAR]", - "The program can also be used with continuous variables (by setting", - "[TT]-nodiscrete[tt]). In this case kinetics of other processes can be", - "studied. This is very much a work in progress and hence the manual", - "(this information) is lagging behind somewhat.[PAR]", - "In order to run [THISMODULE], GROMACS must be compiled with the GNU", - "scientific library." - }; - static int nreplica = 1; - static real tref = 298.15; - static real cutoff = 0.2; - static real ucut = 0.0; - static real Euf = 10; - static real Efu = 30; - static real Ei = 10; - static gmx_bool bHaveT = TRUE; - static real t0 = -1; - static real t1 = -1; - static real tb = 0; - static real te = 0; - static real tol = 1e-3; - static int maxiter = 100; - static int skip = 0; - static int nmult = 1; - static gmx_bool bBack = TRUE; - static gmx_bool bSplit = TRUE; - static gmx_bool bSum = TRUE; - static gmx_bool bDiscrete = TRUE; - t_pargs pa[] = { - { "-time", FALSE, etBOOL, {&bHaveT}, - "Expect a time in the input" }, - { "-b", FALSE, etREAL, {&tb}, - "First time to read from set" }, - { "-e", FALSE, etREAL, {&te}, - "Last time to read from set" }, - { "-bfit", FALSE, etREAL, {&t0}, - "Time to start the fit from" }, - { "-efit", FALSE, etREAL, {&t1}, - "Time to end the fit" }, - { "-T", FALSE, etREAL, {&tref}, - "Reference temperature for computing rate constants" }, - { "-n", FALSE, etINT, {&nreplica}, - "Read data for this number of replicas. Only necessary when files are written in xmgrace format using @type and & as delimiters." }, - { "-cut", FALSE, etREAL, {&cutoff}, - "Cut-off (max) value for regarding a structure as folded" }, - { "-ucut", FALSE, etREAL, {&ucut}, - "Cut-off (max) value for regarding a structure as intermediate (if not folded)" }, - { "-euf", FALSE, etREAL, {&Euf}, - "Initial guess for energy of activation for folding (kJ/mol)" }, - { "-efu", FALSE, etREAL, {&Efu}, - "Initial guess for energy of activation for unfolding (kJ/mol)" }, - { "-ei", FALSE, etREAL, {&Ei}, - "Initial guess for energy of activation for intermediates (kJ/mol)" }, - { "-maxiter", FALSE, etINT, {&maxiter}, - "Max number of iterations" }, - { "-back", FALSE, etBOOL, {&bBack}, - "Take the back reaction into account" }, - { "-tol", FALSE, etREAL, {&tol}, - "Absolute tolerance for convergence of the Nelder and Mead simplex algorithm" }, - { "-skip", FALSE, etINT, {&skip}, - "Skip points in the output [TT].xvg[tt] file" }, - { "-split", FALSE, etBOOL, {&bSplit}, - "Estimate error by splitting the number of replicas in two and refitting" }, - { "-sum", FALSE, etBOOL, {&bSum}, - "Average folding before computing [GRK]chi[grk]^2" }, - { "-discrete", FALSE, etBOOL, {&bDiscrete}, - "Use a discrete folding criterion (F <-> U) or a continuous one" }, - { "-mult", FALSE, etINT, {&nmult}, - "Factor to multiply the data with before discretization" } - }; -#define NPA asize(pa) - - FILE *fp; - real dt_t, dt_d, dt_d2; - int nset_t, nset_d, nset_d2, n_t, n_d, n_d2, i; - const char *tfile, *dfile, *dfile2; - t_remd_data remd; - output_env_t oenv; - - t_filenm fnm[] = { - { efXVG, "-f", "temp", ffREAD }, - { efXVG, "-d", "data", ffREAD }, - { efXVG, "-d2", "data2", ffOPTRD }, - { efXVG, "-o", "ft_all", ffWRITE }, - { efXVG, "-o2", "it_all", ffOPTWR }, - { efXVG, "-o3", "ft_repl", ffOPTWR }, - { efXVG, "-ee", "err_est", ffOPTWR }, - { efLOG, "-g", "remd", ffWRITE }, - { efXVG, "-m", "melt", ffWRITE } - }; -#define NFILE asize(fnm) - - if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_BE_NICE | PCA_TIME_UNIT, - NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv)) - { - return 0; - } - -#ifdef HAVE_LIBGSL - please_cite(stdout, "Spoel2006d"); - if (cutoff < 0) - { - gmx_fatal(FARGS, "cutoff should be >= 0 (rather than %f)", cutoff); - } - - tfile = opt2fn("-f", NFILE, fnm); - dfile = opt2fn("-d", NFILE, fnm); - dfile2 = opt2fn_null("-d2", NFILE, fnm); - - fp = gmx_ffopen(opt2fn("-g", NFILE, fnm), "w"); - - remd.temp = read_xvg_time(tfile, bHaveT, - opt2parg_bSet("-b", NPA, pa), tb, - opt2parg_bSet("-e", NPA, pa), te, - nreplica, &nset_t, &n_t, &dt_t, &remd.time); - printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_t, n_t, tfile, dt_t); - sfree(remd.time); - - remd.data = read_xvg_time(dfile, bHaveT, - opt2parg_bSet("-b", NPA, pa), tb, - opt2parg_bSet("-e", NPA, pa), te, - nreplica, &nset_d, &n_d, &dt_d, &remd.time); - printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_d, n_d, dfile, dt_d); - - if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d)) - { - gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file", - tfile, dfile); - } - - if (dfile2) - { - remd.data2 = read_xvg_time(dfile2, bHaveT, - opt2parg_bSet("-b", NPA, pa), tb, - opt2parg_bSet("-e", NPA, pa), te, - nreplica, &nset_d2, &n_d2, &dt_d2, &remd.time); - printf("Read %d sets of %d points in %s, dt = %g\n\n", - nset_d2, n_d2, dfile2, dt_d2); - if ((nset_d2 != nset_d) || (n_d != n_d2) || (dt_d != dt_d2)) - { - gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file", - dfile, dfile2); - } - } - else - { - remd.data2 = NULL; - } - - remd.nreplica = nset_d; - remd.nframe = n_d; - remd.dt = 1; - preprocess_remd(fp, &remd, cutoff, tref, ucut, bBack, Euf, Efu, Ei, t0, t1, - bSum, bDiscrete, nmult); - - optimize_remd_parameters(&remd, maxiter, tol); - - dump_remd_parameters(fp, &remd, opt2fn("-o", NFILE, fnm), - opt2fn_null("-o2", NFILE, fnm), - opt2fn_null("-o3", NFILE, fnm), - opt2fn_null("-ee", NFILE, fnm), - opt2fn("-m", NFILE, fnm), skip, tref, oenv); - - if (bSplit) - { - printf("Splitting set of replicas in two halves\n"); - for (i = 0; (i < remd.nreplica); i++) - { - remd.bMask[i] = FALSE; - } - remd.nmask = 0; - for (i = 0; (i < remd.nreplica); i += 2) - { - remd.bMask[i] = TRUE; - remd.nmask++; - } - sum_ft(&remd); - optimize_remd_parameters(&remd, maxiter, tol); - dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv); - - for (i = 0; (i < remd.nreplica); i++) - { - remd.bMask[i] = !remd.bMask[i]; - } - remd.nmask = remd.nreplica - remd.nmask; - - sum_ft(&remd); - optimize_remd_parameters(&remd, maxiter, tol); - dump_remd_parameters(fp, &remd, "test2.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv); - - for (i = 0; (i < remd.nreplica); i++) - { - remd.bMask[i] = FALSE; - } - remd.nmask = 0; - for (i = 0; (i < remd.nreplica/2); i++) - { - remd.bMask[i] = TRUE; - remd.nmask++; - } - sum_ft(&remd); - optimize_remd_parameters(&remd, maxiter, tol); - dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv); - - for (i = 0; (i < remd.nreplica); i++) - { - remd.bMask[i] = FALSE; - } - remd.nmask = 0; - for (i = remd.nreplica/2; (i < remd.nreplica); i++) - { - remd.bMask[i] = TRUE; - remd.nmask++; - } - sum_ft(&remd); - optimize_remd_parameters(&remd, maxiter, tol); - dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv); - } - gmx_ffclose(fp); - - view_all(oenv, NFILE, fnm); - -#else - fprintf(stderr, "This program should be compiled with the GNU scientific library. Please install the library and reinstall GROMACS.\n"); -#endif /*HAVE_LIBGSL*/ - - return 0; -} diff --git a/src/gromacs/libgromacs.pc.cmakein b/src/gromacs/libgromacs.pc.cmakein index 082e07f8db..c88b1cb20a 100644 --- a/src/gromacs/libgromacs.pc.cmakein +++ b/src/gromacs/libgromacs.pc.cmakein @@ -5,7 +5,7 @@ Name: libgromacs@GMX_LIBS_SUFFIX@ Description: Gromacs library URL: http://www.gromacs.org Version: @PROJECT_VERSION@ -Requires: @PKG_FFT@ @PKG_XML@ @PKG_GSL@ +Requires: @PKG_FFT@ @PKG_XML@ Libs.private: @CMAKE_THREAD_LIBS_INIT@ @PKG_DL_LIBS@ @OpenMP_LINKER_FLAGS@ Libs: -L${libdir} -lgromacs@GMX_LIBS_SUFFIX@ @PKG_FFT_LIBS@ -lm Cflags: -I${includedir} @PKG_CFLAGS@ diff --git a/src/programs/CreateLinks.cmake.cmakein b/src/programs/CreateLinks.cmake.cmakein index 5e663d2ff7..04a4905438 100644 --- a/src/programs/CreateLinks.cmake.cmakein +++ b/src/programs/CreateLinks.cmake.cmakein @@ -90,7 +90,6 @@ set(SYMLINK_NAMES g_helix g_helixorient g_hydorder - g_kinetics g_lie g_mdmat g_mindist diff --git a/src/programs/legacymodules.cpp b/src/programs/legacymodules.cpp index e7ba8b86df..8b014b2663 100644 --- a/src/programs/legacymodules.cpp +++ b/src/programs/legacymodules.cpp @@ -261,8 +261,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager) "Calculate local pitch/bending/rotation/orientation inside helices"); registerModule(manager, &gmx_hydorder, "hydorder", "Compute tetrahedrality parameters around a given atom"); - registerModule(manager, &gmx_kinetics, "kinetics", - "Analyze kinetic constants from properties based on the Eyring model"); registerModule(manager, &gmx_lie, "lie", "Estimate free energy from linear combinations"); registerModule(manager, &gmx_mdmat, "mdmat", @@ -459,7 +457,6 @@ void registerLegacyModules(gmx::CommandLineModuleManager *manager) group.addModule("current"); group.addModule("dos"); group.addModule("dyecoupl"); - group.addModule("kinetics"); group.addModule("principal"); group.addModule("tcaf"); group.addModule("traj"); -- 2.22.0