add_definitions(-DTMPI_ATOMICS_DISABLED)
endif()
- # now that we have detected the dependencies, do the second GPU configure pass
- gmx_gpu_setup()
+# Note this relies on zlib detection having already run
+include(gmxManageTNG)
+
+ if(GMX_GPU)
+ # now that we have detected the dependencies, do the second configure pass
+ gmx_gpu_setup()
+ else()
+ mark_as_advanced(CUDA_HOST_COMPILER)
+ endif()
if(CYGWIN)
set(GMX_CYGWIN 1)
message(STATUS "Detecting best SIMD instructions for this CPU")
# Get CPU SIMD properties information
- set(_compile_definitions "${GCC_INLINE_ASM_DEFINE} -I${CMAKE_SOURCE_DIR}/src -DGMX_CPUID_STANDALONE")
+ set(_compile_definitions "${GCC_INLINE_ASM_DEFINE} -I${CMAKE_SOURCE_DIR}/src -DGMX_CPUINFO_STANDALONE ${GMX_STDLIB_CXX_FLAGS}")
+ if(GMX_TARGET_X86)
+ set(_compile_definitions "${_compile_definitions} -DGMX_TARGET_X86")
+ endif()
+
+ # Prepare a default suggestion
+ set(OUTPUT_SIMD "None")
# We need to execute the binary, so this only works if not cross-compiling.
# However, note that we are NOT limited to x86.
if(NOT CMAKE_CROSSCOMPILING)
- # We can use try_run(... LINK_LIBRARIES ${GMX_STDLIB_LIBRARIES})
- # once we require CMake at least 2.8.11. The code below works
- # even when GMX_STDLIB_LIBRARIES is empty.
+ # TODO Extract this try_compile to a helper function, because
+ # it duplicates code in gmxDetectSimd.cmake
+ set(GMX_DETECTSIMD_BINARY "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/GmxDetectSimd${CMAKE_EXECUTABLE_SUFFIX}")
+ set(LINK_LIBRARIES "${GMX_STDLIB_LIBRARIES}")
- try_run(GMX_CPUINFO_RUN_SIMD GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- CMAKE_FLAGS "-DLINK_LIBRARIES=${LINK_LIBRARIES}"
- RUN_OUTPUT_VARIABLE OUTPUT_TMP
- COMPILE_OUTPUT_VARIABLE GMX_CPUINFO_COMPILE_OUTPUT
- ARGS "-features")
-
- if(NOT GMX_CPUINFO_COMPILED)
- message(WARNING "Cannot compile cpuinfo code, which means no SIMD instructions.")
- message(STATUS "Compile output: ${GMX_CPUINFO_COMPILE_OUTPUT}")
- set(OUTPUT_TMP "None")
- elseif(NOT GMX_CPUINFO_RUN_SIMD EQUAL 0)
- message(WARNING "Cannot run cpuinfo code, which means no SIMD instructions.")
- message(STATUS "Run output: ${OUTPUT_TMP}")
- set(OUTPUT_TMP "None")
- endif(NOT GMX_CPUINFO_COMPILED)
+ try_compile(GMX_DETECTSIMD_COMPILED
+ "${CMAKE_CURRENT_BINARY_DIR}"
- "${CMAKE_CURRENT_SOURCE_DIR}/src/gromacs/gmxlib/gmx_cpuid.c"
++ "${CMAKE_CURRENT_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp"
+ COMPILE_DEFINITIONS "${_compile_definitions}"
++ CMAKE_FLAGS "-DLINK_LIBRARIES=${LINK_LIBRARIES}"
+ OUTPUT_VARIABLE GMX_DETECTSIMD_COMPILED_OUTPUT
+ COPY_FILE ${GMX_DETECTSIMD_BINARY})
+ unset(_compile_definitions)
- set(OUTPUT_SIMD "None")
- if(GMX_TARGET_X86)
- if(OUTPUT_TMP MATCHES " avx512er ")
- set(OUTPUT_SIMD "AVX_512_KNL")
- elseif(OUTPUT_TMP MATCHES " avx512f ")
- set(OUTPUT_SIMD "AVX_512")
- elseif(OUTPUT_TMP MATCHES " avx2 ")
- set(OUTPUT_SIMD "AVX2_256")
- elseif(OUTPUT_TMP MATCHES " avx ")
- if(OUTPUT_TMP MATCHES " fma4 ")
- # AMD that works better with avx-128-fma
- set(OUTPUT_SIMD "AVX_128_FMA")
+ if(GMX_DETECTSIMD_COMPILED)
+ # TODO Extract this duplication of
+ # gmxSetBuildInformation.cmake to a helper function
+ if(NOT DEFINED GMX_DETECTSIMD_RUN)
+ execute_process(COMMAND ${GMX_DETECTSIMD_BINARY} "-simd"
+ RESULT_VARIABLE GMX_DETECTSIMD_RUN
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_DETECTSIMD_RUN "${GMX_DETECTSIMD_RUN}" CACHE INTERNAL "Result of running CPUID code with arg -simd")
+ if(GMX_DETECTSIMD_RUN EQUAL 0)
+ # Make a concrete suggestion of SIMD level
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_SIMD)
+ message(STATUS "Detected best SIMD instructions for this CPU - ${OUTPUT_SIMD}")
else()
- # Intel
- set(OUTPUT_SIMD "AVX_256")
+ message(WARNING "Cannot run CPUID code, which means no SIMD suggestion can be made.")
+ message(STATUS "Run output: ${OUTPUT_TMP}")
endif()
- elseif(OUTPUT_TMP MATCHES " sse4.1 ")
- set(OUTPUT_SIMD "SSE4.1")
- elseif(OUTPUT_TMP MATCHES " sse2 ")
- set(OUTPUT_SIMD "SSE2")
endif()
else()
- if(OUTPUT_TMP MATCHES " vsx ")
- set(OUTPUT_SIMD "IBM_VSX")
- elseif(OUTPUT_TMP MATCHES " vmx ")
- set(OUTPUT_SIMD "IBM_VMX")
- elseif(OUTPUT_TMP MATCHES " qpx ")
- set(OUTPUT_SIMD "IBM_QPX")
- elseif(OUTPUT_TMP MATCHES " neon_asimd ")
- set(OUTPUT_SIMD "ARM_NEON_ASIMD")
- elseif(OUTPUT_TMP MATCHES " neon ")
- set(OUTPUT_SIMD "ARM_NEON")
- endif()
+ message(WARNING "Cannot compile CPUID code, which means no SIMD instructions.")
+ message(STATUS "Compile output: ${GMX_CPUID_COMPILE_OUTPUT}")
endif()
-
- set(${_suggested_simd} "${OUTPUT_SIMD}" PARENT_SCOPE)
- message(STATUS "Detected best SIMD instructions for this CPU - ${OUTPUT_SIMD}")
else()
- set(${_suggested_simd} "None" PARENT_SCOPE)
message(WARNING "Cannot detect SIMD architecture for this cross-compile; you should check it manually.")
endif()
+
+ set(${_suggested_simd} "${OUTPUT_SIMD}" CACHE INTERNAL "Suggested SIMD")
endfunction()
function(gmx_detect_simd _suggested_simd)
message(STATUS "Setting build user & time - not on Unix, using anonymous")
endif()
- if(NOT CMAKE_CROSSCOMPILING)
- # Get CPU information, e.g. for deciding what SIMD support exists
- set(_compile_definitions "${GCC_INLINE_ASM_DEFINE} -I${CMAKE_SOURCE_DIR}/src -DGMX_CPUINFO_STANDALONE")
- try_run(GMX_CPUINFO_RUN_VENDOR GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_VENDOR ARGS "-vendor")
- try_run(GMX_CPUINFO_RUN_BRAND GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_BRAND ARGS "-brand")
- try_run(GMX_CPUINFO_RUN_FAMILY GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_FAMILY ARGS "-family")
- try_run(GMX_CPUINFO_RUN_MODEL GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_MODEL ARGS "-model")
- try_run(GMX_CPUINFO_RUN_STEPPING GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_STEPPING ARGS "-stepping")
- try_run(GMX_CPUINFO_RUN_FEATURES GMX_CPUINFO_COMPILED
- ${CMAKE_BINARY_DIR}
- ${CMAKE_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp
- COMPILE_DEFINITIONS ${_compile_definitions}
- RUN_OUTPUT_VARIABLE OUTPUT_CPU_FEATURES ARGS "-features")
- unset(_compile_definitions)
-
- string(STRIP "${OUTPUT_CPU_VENDOR}" OUTPUT_CPU_VENDOR)
- string(STRIP "${OUTPUT_CPU_BRAND}" OUTPUT_CPU_BRAND)
- string(STRIP "${OUTPUT_CPU_FAMILY}" OUTPUT_CPU_FAMILY)
- string(STRIP "${OUTPUT_CPU_MODEL}" OUTPUT_CPU_MODEL)
- string(STRIP "${OUTPUT_CPU_STEPPING}" OUTPUT_CPU_STEPPING)
- string(STRIP "${OUTPUT_CPU_FEATURES}" OUTPUT_CPU_FEATURES)
+ # Set up some defaults that will usually be overridden
+ if(CMAKE_CROSSCOMPILING)
+ set(_reason ", cross-compiled")
+ endif()
+ set(OUTPUT_CPU_VENDOR "Unknown${_reason}")
+ set(OUTPUT_CPU_BRAND "Unknown${_reason}")
+ set(OUTPUT_CPU_FAMILY "0")
+ set(OUTPUT_CPU_MODEL "0")
+ set(OUTPUT_CPU_STEPPING "0")
+ set(OUTPUT_CPU_FEATURES "Unknown${_reason}")
+ unset(_reason)
- if(GMX_CPUINFO_RUN_VENDOR EQUAL 0)
- set(BUILD_CPU_VENDOR "${OUTPUT_CPU_VENDOR}" CACHE INTERNAL "Build CPU vendor")
- else()
- set(BUILD_CPU_VENDOR "Unknown, detect failed" CACHE INTERNAL "Build CPU vendor")
- endif()
- if(GMX_CPUINFO_RUN_BRAND EQUAL 0)
- set(BUILD_CPU_BRAND "${OUTPUT_CPU_BRAND}" CACHE INTERNAL "Build CPU brand")
- else()
- set(BUILD_CPU_BRAND "Unknown, detect failed" CACHE INTERNAL "Build CPU brand")
- endif()
- if(GMX_CPUINFO_RUN_FAMILY EQUAL 0)
- set(BUILD_CPU_FAMILY "${OUTPUT_CPU_FAMILY}" CACHE INTERNAL "Build CPU family")
- else()
- set(BUILD_CPU_FAMILY "0" CACHE INTERNAL "Build CPU family")
- endif()
- if(GMX_CPUINFO_RUN_MODEL EQUAL 0)
- set(BUILD_CPU_MODEL "${OUTPUT_CPU_MODEL}" CACHE INTERNAL "Build CPU model")
- else()
- set(BUILD_CPU_MODEL "0" CACHE INTERNAL "Build CPU model")
- endif()
- if(GMX_CPUINFO_RUN_STEPPING EQUAL 0)
- set(BUILD_CPU_STEPPING "${OUTPUT_CPU_STEPPING}" CACHE INTERNAL "Build CPU stepping")
- else()
- set(BUILD_CPU_STEPPING "0" CACHE INTERNAL "Build CPU stepping")
- endif()
- if(GMX_CPUINFO_RUN_FEATURES EQUAL 0)
- set(BUILD_CPU_FEATURES "${OUTPUT_CPU_FEATURES}" CACHE INTERNAL "Build CPU features")
- else()
- set(BUILD_CPU_FEATURES "" CACHE INTERNAL "Build CPU features")
+ if(NOT CMAKE_CROSSCOMPILING)
+ # Get CPU information, e.g. for deciding what SIMD support probably exists
+ set(_compile_definitions "${GCC_INLINE_ASM_DEFINE} -I${CMAKE_SOURCE_DIR}/src -DGMX_CPUID_STANDALONE")
+ if(GMX_TARGET_X86)
+ set(_compile_definitions "${_compile_definitions} -DGMX_TARGET_X86")
endif()
- else()
-
- set(BUILD_CPU_VENDOR "Unknown, cross-compiled" CACHE INTERNAL "Build CPU vendor")
- set(BUILD_CPU_BRAND "Unknown, cross-compiled" CACHE INTERNAL "Build CPU brand")
- set(BUILD_CPU_FAMILY "0" CACHE INTERNAL "Build CPU family")
- set(BUILD_CPU_MODEL "0" CACHE INTERNAL "Build CPU model")
- set(BUILD_CPU_STEPPING "0" CACHE INTERNAL "Build CPU stepping")
- set(BUILD_CPU_FEATURES "" CACHE INTERNAL "Build CPU features")
+ set(GMX_BUILDINFORMATION_BINARY "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/GmxBuildInformation${CMAKE_EXECUTABLE_SUFFIX}")
+ # TODO Extract this try_compile to a helper function, because
+ # it duplicates code in gmxDetectSimd.cmake
+ try_compile(GMX_BUILDINFORMATION_COMPILED
+ "${CMAKE_CURRENT_BINARY_DIR}"
- "${CMAKE_CURRENT_SOURCE_DIR}/src/gromacs/gmxlib/gmx_cpuid.c"
++ "${CMAKE_CURRENT_SOURCE_DIR}/src/gromacs/hardware/cpuinfo.cpp"
+ COMPILE_DEFINITIONS "${_compile_definitions}"
+ OUTPUT_VARIABLE GMX_BUILDINFORMATION_COMPILED_OUTPUT
+ COPY_FILE ${GMX_BUILDINFORMATION_BINARY})
+ unset(_compile_definitions)
+ if(GMX_BUILDINFORMATION_COMPILED)
+ # TODO Extract this duplication to a helper function (also
+ # from gmxDetectSimd.cmake)
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_VENDOR)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-vendor"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_VENDOR
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_VENDOR "${GMX_BUILDINFORMATION_RUN_VENDOR}" CACHE INTERNAL "Result of running CPUID code with arg -vendor")
+ if(GMX_BUILDINFORMATION_RUN_VENDOR EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_VENDOR)
+ endif()
+ endif()
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_BRAND)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-brand"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_BRAND
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_BRAND "${GMX_BUILDINFORMATION_RUN_BRAND}" CACHE INTERNAL "Result of running CPUID code with arg -brand")
+ if(GMX_BUILDINFORMATION_RUN_BRAND EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_BRAND)
+ endif()
+ endif()
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_FAMILY)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-family"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_FAMILY
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_FAMILY "${GMX_BUILDINFORMATION_RUN_FAMILY}" CACHE INTERNAL "Result of running CPUID code with arg -family")
+ if(GMX_BUILDINFORMATION_RUN_FAMILY EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_FAMILY)
+ endif()
+ endif()
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_MODEL)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-model"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_MODEL
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_MODEL "${GMX_BUILDINFORMATION_RUN_MODEL}" CACHE INTERNAL "Result of running CPUID code with arg -model")
+ if(GMX_BUILDINFORMATION_RUN_MODEL EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_MODEL)
+ endif()
+ endif()
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_STEPPING)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-stepping"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_STEPPING
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_STEPPING "${GMX_BUILDINFORMATION_RUN_STEPPING}" CACHE INTERNAL "Result of running CPUID code with arg -stepping")
+ if(GMX_BUILDINFORMATION_RUN_STEPPING EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_STEPPING)
+ endif()
+ endif()
+ if(NOT DEFINED GMX_BUILDINFORMATION_RUN_FEATURES)
+ execute_process(COMMAND ${GMX_BUILDINFORMATION_BINARY} "-features"
+ RESULT_VARIABLE GMX_BUILDINFORMATION_RUN_FEATURES
+ OUTPUT_VARIABLE OUTPUT_TMP
+ ERROR_QUIET)
+ set(GMX_BUILDINFORMATION_RUN_FEATURES "${GMX_BUILDINFORMATION_RUN_FEATURES}" CACHE INTERNAL "Result of running CPUID code with arg -features")
+ if(GMX_BUILDINFORMATION_RUN_FEATURES EQUAL 0)
+ string(STRIP "${OUTPUT_TMP}" OUTPUT_CPU_FEATURES)
+ endif()
+ endif()
+ endif()
endif()
+ set(BUILD_CPU_VENDOR "${OUTPUT_CPU_VENDOR}" CACHE INTERNAL "Build CPU vendor")
+ set(BUILD_CPU_BRAND "${OUTPUT_CPU_BRAND}" CACHE INTERNAL "Build CPU brand")
+ set(BUILD_CPU_FAMILY "${OUTPUT_CPU_FAMILY}" CACHE INTERNAL "Build CPU family")
+ set(BUILD_CPU_MODEL "${OUTPUT_CPU_MODEL}" CACHE INTERNAL "Build CPU model")
+ set(BUILD_CPU_STEPPING "${OUTPUT_CPU_STEPPING}" CACHE INTERNAL "Build CPU stepping")
+ set(BUILD_CPU_FEATURES "${OUTPUT_CPU_FEATURES}" CACHE INTERNAL "Build CPU features")
+
ENDIF(NOT DEFINED BUILD_USER)
endmacro(gmx_set_build_information)
--- /dev/null
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+/*
+ * 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,2015,2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "filenm.h"
+
+#include <cstdio>
+#include <cstring>
+
+#include "gromacs/fileio/filetypes.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/smalloc.h"
+
+/* Use bitflag ... */
+#define IS_SET(fn) ((fn.flag &ffSET) != 0)
+#define IS_OPT(fn) ((fn.flag &ffOPT) != 0)
+
++const t_filenm *getFilenm(const char *opt, int nfile, const t_filenm fnm[])
++{
++ int i;
++
++ for (i = 0; (i < nfile); i++)
++ {
++ if (strcmp(opt, fnm[i].opt) == 0)
++ {
++ return &fnm[i];
++ }
++ }
++
++ return NULL;
++}
++
+const char *opt2fn(const char *opt, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (std::strcmp(opt, fnm[i].opt) == 0)
+ {
+ return fnm[i].fns[0];
+ }
+ }
+
+ fprintf(stderr, "No option %s\n", opt);
+
+ return NULL;
+}
+
+int opt2fns(char **fns[], const char *opt, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (strcmp(opt, fnm[i].opt) == 0)
+ {
+ *fns = fnm[i].fns;
+ return fnm[i].nfiles;
+ }
+ }
+
+ fprintf(stderr, "No option %s\n", opt);
+ return 0;
+}
+
+const char *ftp2fn(int ftp, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (ftp == fnm[i].ftp)
+ {
+ return fnm[i].fns[0];
+ }
+ }
+
+ fprintf(stderr, "ftp2fn: No filetype %s\n", ftp2ext_with_dot(ftp));
+ return NULL;
+}
+
+int ftp2fns(char **fns[], int ftp, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (ftp == fnm[i].ftp)
+ {
+ *fns = fnm[i].fns;
+ return fnm[i].nfiles;
+ }
+ }
+
+ fprintf(stderr, "ftp2fn: No filetype %s\n", ftp2ext_with_dot(ftp));
+ return 0;
+}
+
+gmx_bool ftp2bSet(int ftp, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (ftp == fnm[i].ftp)
+ {
+ return (gmx_bool) IS_SET(fnm[i]);
+ }
+ }
+
+ fprintf(stderr, "ftp2fn: No filetype %s\n", ftp2ext_with_dot(ftp));
+
+ return FALSE;
+}
+
+gmx_bool opt2bSet(const char *opt, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (std::strcmp(opt, fnm[i].opt) == 0)
+ {
+ return (gmx_bool) IS_SET(fnm[i]);
+ }
+ }
+
+ fprintf(stderr, "No option %s\n", opt);
+
+ return FALSE;
+}
+
+const char *opt2fn_null(const char *opt, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (std::strcmp(opt, fnm[i].opt) == 0)
+ {
+ if (IS_OPT(fnm[i]) && !IS_SET(fnm[i]))
+ {
+ return NULL;
+ }
+ else
+ {
+ return fnm[i].fns[0];
+ }
+ }
+ }
+ fprintf(stderr, "No option %s\n", opt);
+ return NULL;
+}
+
+const char *ftp2fn_null(int ftp, int nfile, const t_filenm fnm[])
+{
+ int i;
+
+ for (i = 0; (i < nfile); i++)
+ {
+ if (ftp == fnm[i].ftp)
+ {
+ if (IS_OPT(fnm[i]) && !IS_SET(fnm[i]))
+ {
+ return NULL;
+ }
+ else
+ {
+ return fnm[i].fns[0];
+ }
+ }
+ }
+ fprintf(stderr, "ftp2fn: No filetype %s\n", ftp2ext_with_dot(ftp));
+ return NULL;
+}
+
+gmx_bool is_optional(const t_filenm *fnm)
+{
+ return ((fnm->flag & ffOPT) == ffOPT);
+}
+
+gmx_bool is_output(const t_filenm *fnm)
+{
+ return ((fnm->flag & ffWRITE) == ffWRITE);
+}
+
+gmx_bool is_set(const t_filenm *fnm)
+{
+ return ((fnm->flag & ffSET) == ffSET);
+}
+
+int add_suffix_to_output_names(t_filenm *fnm, int nfile, const char *suffix)
+{
+ int i, j;
+ char buf[STRLEN], newname[STRLEN];
+ char *extpos;
+
+ for (i = 0; i < nfile; i++)
+ {
+ if (is_output(&fnm[i]) && fnm[i].ftp != efCPT)
+ {
+ /* We never use multiple _outputs_, but we might as well check
+ for it, just in case... */
+ for (j = 0; j < fnm[i].nfiles; j++)
+ {
+ std::strncpy(buf, fnm[i].fns[j], STRLEN - 1);
+ extpos = strrchr(buf, '.');
+ *extpos = '\0';
+ sprintf(newname, "%s%s.%s", buf, suffix, extpos + 1);
+ sfree(fnm[i].fns[j]);
+ fnm[i].fns[j] = gmx_strdup(newname);
+ }
+ }
+ }
+ return 0;
+}
+
+t_filenm *dup_tfn(int nf, const t_filenm tfn[])
+{
+ int i, j;
+ t_filenm *ret;
+
+ snew(ret, nf);
+ for (i = 0; i < nf; i++)
+ {
+ ret[i] = tfn[i]; /* just directly copy all non-string fields */
+ if (tfn[i].opt)
+ {
+ ret[i].opt = gmx_strdup(tfn[i].opt);
+ }
+ else
+ {
+ ret[i].opt = NULL;
+ }
+
+ if (tfn[i].fn)
+ {
+ ret[i].fn = gmx_strdup(tfn[i].fn);
+ }
+ else
+ {
+ ret[i].fn = NULL;
+ }
+
+ if (tfn[i].nfiles > 0)
+ {
+ snew(ret[i].fns, tfn[i].nfiles);
+ for (j = 0; j < tfn[i].nfiles; j++)
+ {
+ ret[i].fns[j] = gmx_strdup(tfn[i].fns[j]);
+ }
+ }
+ }
+ return ret;
+}
+
+void done_filenms(int nf, t_filenm fnm[])
+{
+ int i, j;
+
+ for (i = 0; i < nf; ++i)
+ {
+ for (j = 0; j < fnm[i].nfiles; ++j)
+ {
+ sfree(fnm[i].fns[j]);
+ }
+ sfree(fnm[i].fns);
+ fnm[i].fns = NULL;
+ }
+}
--- /dev/null
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+/*
+ * 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,2015,2016, 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.
+ */
+/*! \file
+ * \brief
+ * Declares t_filenm for old-style command-line parsing of file name options.
+ *
+ * \inpublicapi
+ * \ingroup module_commandline
+ */
+#ifndef GMX_COMMANDLINE_FILENM_H
+#define GMX_COMMANDLINE_FILENM_H
+
+#include "gromacs/fileio/filetypes.h"
+#include "gromacs/utility/basedefinitions.h"
+
+struct t_commrec;
+
+//! \addtogroup module_commandline
+//! \{
+
+/*! \brief
+ * File name option definition for C code.
+ *
+ * \inpublicapi
+ */
+struct t_filenm {
+ int ftp; //!< File type (see enum in filetypes.h)
+ const char *opt; //!< Command line option
+ const char *fn; //!< File name (as set in source code)
+ unsigned long flag; //!< Flag for all kinds of info (see defs)
+ int nfiles; //!< number of files
+ char **fns; //!< File names
+};
+
+//! Whether a file name option is set.
+#define ffSET 1<<0
+//! Whether a file name option specifies an input file.
+#define ffREAD 1<<1
+//! Whether a file name option specifies an output file.
+#define ffWRITE 1<<2
+//! Whether a file name option specifies an optional file.
+#define ffOPT 1<<3
+//! Whether a file name option specifies a library file.
+#define ffLIB 1<<4
+//! Whether a file name option accepts multiple file names.
+#define ffMULT 1<<5
+//! Whether an input file name option accepts non-existent files.
+#define ffALLOW_MISSING 1<<6
+//! Convenience flag for an input/output file.
+#define ffRW (ffREAD | ffWRITE)
+//! Convenience flag for an optional input file.
+#define ffOPTRD (ffREAD | ffOPT)
+//! Convenience flag for an optional output file.
+#define ffOPTWR (ffWRITE| ffOPT)
+//! Convenience flag for an optional input/output file.
+#define ffOPTRW (ffRW | ffOPT)
+//! Convenience flag for a library input file.
+#define ffLIBRD (ffREAD | ffLIB)
+//! Convenience flag for an optional library input file.
+#define ffLIBOPTRD (ffOPTRD | ffLIB)
+//! Convenience flag for an input file that accepts multiple files.
+#define ffRDMULT (ffREAD | ffMULT)
+//! Convenience flag for an optional input file that accepts multiple files.
+#define ffOPTRDMULT (ffRDMULT | ffOPT)
+//! Convenience flag for an output file that accepts multiple files.
+#define ffWRMULT (ffWRITE | ffMULT)
+//! Convenience flag for an optional output file that accepts multiple files.
+#define ffOPTWRMULT (ffWRMULT | ffOPT)
+
+/*! \brief
+ * Returns the filename belonging to cmd-line option opt, or NULL when
+ * no such option.
+ */
+const char *opt2fn(const char *opt, int nfile, const t_filenm fnm[]);
+
+/*! \brief
+ * Returns the filenames belonging to cmd-line option opt, or NULL when
+ * no such option.
+ */
+int opt2fns(char **fns[], const char *opt, int nfile,
+ const t_filenm fnm[]);
+
++/*! \brief
++ * Return a pointer to the t_filenm data structure of filenames belonging to
++ * command-line option opt, or NULL when no such option was used.
++ */
++const t_filenm *getFilenm(const char *opt, int nfile, const t_filenm fnm[]);
++
+//! Returns a file pointer from the filename.
+#define opt2FILE(opt, nfile, fnm, mode) gmx_ffopen(opt2fn(opt, nfile, fnm), mode)
+
+//! Returns the first file name with type ftp, or NULL when none found.
+const char *ftp2fn(int ftp, int nfile, const t_filenm fnm[]);
+
+/*! \brief
+ * Returns the number of files for the first option with type ftp
+ * and the files in **fns[] (will be allocated), or NULL when none found.
+ */
+int ftp2fns(char **fns[], int ftp, int nfile, const t_filenm fnm[]);
+
+//! Returns a file pointer from the file type.
+#define ftp2FILE(ftp, nfile, fnm, mode) gmx_ffopen(ftp2fn(ftp, nfile, fnm), mode)
+
+//! Returns TRUE when this file type has been found on the cmd-line.
+gmx_bool ftp2bSet(int ftp, int nfile, const t_filenm fnm[]);
+
+//! Returns TRUE when this option has been found on the cmd-line.
+gmx_bool opt2bSet(const char *opt, int nfile, const t_filenm fnm[]);
+
+/*! \brief
+ * Returns the file name belonging top cmd-line option opt, or NULL when
+ * no such option.
+ *
+ * Also return NULL when opt is optional and option is not set.
+ */
+const char *opt2fn_null(const char *opt, int nfile, const t_filenm fnm[]);
+
+/*! \brief
+ * Returns the first file name with type ftp, or NULL when none found.
+ *
+ * Also return NULL when ftp is optional and option is not set.
+ */
+const char *ftp2fn_null(int ftp, int nfile, const t_filenm fnm[]);
+
+//! Returns whether or not this filenm is optional.
+gmx_bool is_optional(const t_filenm *fnm);
+
+//! Returns whether or not this filenm is output.
+gmx_bool is_output(const t_filenm *fnm);
+
+//! Returns whether or not this filenm is set.
+gmx_bool is_set(const t_filenm *fnm);
+
+/*! \brief
+ * When we do checkpointing, this routine is called to check for previous
+ * output files and append a '.partNNNN' suffix before the (output) file extensions.
+ */
+int add_suffix_to_output_names(t_filenm *fnm, int nfile, const char *suffix);
+
+/*! \brief
+ * Duplicates the filename list (to make a private copy for each thread,
+ * for example).
+ */
+t_filenm *dup_tfn(int nf, const t_filenm tfn[]);
+
+//! Frees memory allocated for file names by parse_common_args().
+void done_filenms(int nf, t_filenm fnm[]);
+
+//! \}
+
+#endif
--- /dev/null
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+/*
+ * 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,2015,2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "filetypes.h"
+
+#include <cstring>
+
+#include "gromacs/utility/arraysize.h"
+#include "gromacs/utility/cstringutil.h"
+
+enum
+{
+ eftASC, eftXDR, eftTNG, eftGEN, eftNR
+};
+
+/* To support multiple file types with one general (eg TRX) we have
+ * these arrays.
+ */
+static const int trxs[] =
+{
+ efXTC, efTRR, efCPT,
+ efGRO, efG96, efPDB, efTNG
+};
+#define NTRXS asize(trxs)
+
+static const int trcompressed[] =
+{
+ efXTC,
+ efTNG
+};
+#define NTRCOMPRESSED asize(trcompressed)
+
+static const int tros[] =
+{
+ efXTC, efTRR,
+ efGRO, efG96, efPDB, efTNG
+};
+#define NTROS asize(tros)
+
+static const int trns[] =
+{
+ efTRR, efCPT,
+ efTNG
+};
+#define NTRNS asize(trns)
+
+static const int stos[] =
+{ efGRO, efG96, efPDB, efBRK, efENT, efESP };
+#define NSTOS asize(stos)
+
+static const int stxs[] =
+{
+ efGRO, efG96, efPDB, efBRK, efENT, efESP,
+ efTPR
+};
+#define NSTXS asize(stxs)
+
+static const int tpss[] =
+{
+ efTPR,
+ efGRO, efG96, efPDB, efBRK, efENT
+};
+#define NTPSS asize(tpss)
+
+typedef struct
+{
+ int ftype;
+ const char *ext;
+ const char *defnm;
+ const char *defopt;
+ const char *descr;
+ int ntps;
+ const int *tps;
+} t_deffile;
+
+/* this array should correspond to the enum in filetypes.h */
+static const t_deffile deffile[efNR] =
+{
+ { eftASC, ".mdp", "grompp", "-f", "grompp input file with MD parameters" },
+ { eftGEN, ".???", "traj", "-f", "Trajectory", NTRXS, trxs },
+ { eftGEN, ".???", "trajout", "-f", "Trajectory", NTROS, tros },
+ { eftGEN, ".???", "traj", NULL,
+ "Full precision trajectory", NTRNS, trns },
+ { eftXDR, ".trr", "traj", NULL, "Trajectory in portable xdr format" },
+ { eftGEN, ".???", "traj_comp", NULL,
+ "Compressed trajectory (tng format or portable xdr format)", NTRCOMPRESSED, trcompressed},
+ { eftXDR, ".xtc", "traj", NULL,
+ "Compressed trajectory (portable xdr format): xtc" },
+ { eftTNG, ".tng", "traj", NULL,
+ "Trajectory file (tng format)" },
+ { eftXDR, ".edr", "ener", NULL, "Energy file"},
+ { eftGEN, ".???", "conf", "-c", "Structure file", NSTXS, stxs },
+ { eftGEN, ".???", "out", "-o", "Structure file", NSTOS, stos },
+ { eftASC, ".gro", "conf", "-c", "Coordinate file in Gromos-87 format" },
+ { eftASC, ".g96", "conf", "-c", "Coordinate file in Gromos-96 format" },
+ { eftASC, ".pdb", "eiwit", "-f", "Protein data bank file"},
+ { eftASC, ".brk", "eiwit", "-f", "Brookhaven data bank file"},
+ { eftASC, ".ent", "eiwit", "-f", "Entry in the protein date bank" },
+ { eftASC, ".esp", "conf", "-f", "Coordinate file in Espresso format" },
+ { eftASC, ".pqr", "state", "-o", "Coordinate file for MEAD"},
+ { eftXDR, ".cpt", "state", "-cp", "Checkpoint file"},
+ { eftASC, ".log", "run", "-l", "Log file"},
+ { eftASC, ".xvg", "graph", "-o", "xvgr/xmgr file"},
+ { eftASC, ".out", "hello", "-o", "Generic output file"},
+ { eftASC, ".ndx", "index", "-n", "Index file", },
+ { eftASC, ".top", "topol", "-p", "Topology file"},
+ { eftASC, ".itp", "topinc", NULL, "Include file for topology"},
+ { eftGEN, ".???", "topol", "-s", "Structure+mass(db)", NTPSS, tpss },
+ { eftXDR, ".tpr", "topol", "-s", "Portable xdr run input file"},
+ { eftASC, ".tex", "doc", "-o", "LaTeX file"},
+ { eftASC, ".rtp", "residue", NULL, "Residue Type file used by pdb2gmx" },
+ { eftASC, ".atp", "atomtp", NULL, "Atomtype file used by pdb2gmx" },
+ { eftASC, ".hdb", "polar", NULL, "Hydrogen data base"},
+ { eftASC, ".dat", "nnnice", NULL, "Generic data file"},
+ { eftASC, ".dlg", "user", NULL, "Dialog Box data for ngmx"},
+ { eftASC, ".map", "ss", NULL, "File that maps matrix data to colors" },
+ { eftASC, ".eps", "plot", NULL, "Encapsulated PostScript (tm) file" },
+ { eftASC, ".mat", "ss", NULL, "Matrix Data file"},
+ { eftASC, ".m2p", "ps", NULL, "Input file for mat2ps"},
+ { eftXDR, ".mtx", "hessian", "-m", "Hessian matrix"},
+ { eftASC, ".edi", "sam", NULL, "ED sampling input"},
+ { eftASC, ".cub", "pot", NULL, "Gaussian cube file" },
+ { eftASC, ".xpm", "root", NULL, "X PixMap compatible matrix file" },
+ { eftASC, "", "rundir", NULL, "Run directory" }
+};
+
+const char *ftp2ext(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].ext[0] != '\0' ? deffile[ftp].ext + 1 : "";
+ }
+ else
+ {
+ return "unknown";
+ }
+}
+
+const char *ftp2ext_generic(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ switch (ftp)
+ {
+ case efTRX:
+ return "trx";
+ case efTRN:
+ return "trn";
+ case efSTO:
+ return "sto";
+ case efSTX:
+ return "stx";
+ case efTPS:
+ return "tps";
+ default:
+ return ftp2ext(ftp);
+ }
+ }
+ else
+ {
+ return "unknown";
+ }
+}
+
+const char *ftp2ext_with_dot(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].ext;
+ }
+ else
+ {
+ return "unknown";
+ }
+}
+
+int ftp2generic_count(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].ntps;
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+const int *ftp2generic_list(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].tps;
+ }
+ else
+ {
+ return 0;
+ }
+}
+
+const char *ftp2desc(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].descr;
+ }
+ else
+ {
+ return "unknown filetype";
+ }
+}
+
+gmx_bool ftp_is_text(int ftp)
+{
+ if ((ftp >= 0) && (ftp < efNR))
+ {
+ return deffile[ftp].ftype == eftASC;
+ }
+ return FALSE;
+}
+
+gmx_bool ftp_is_xdr(int ftp)
+{
+ if ((ftp >= 0) && (ftp < efNR))
+ {
+ return deffile[ftp].ftype == eftXDR;
+ }
+ return FALSE;
+}
+
+const char *ftp2defnm(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].defnm;
+ }
+ else
+ {
+ return NULL;
+ }
+}
+
+const char *ftp2defopt(int ftp)
+{
+ if ((0 <= ftp) && (ftp < efNR))
+ {
+ return deffile[ftp].defopt;
+ }
+ else
+ {
+ return NULL;
+ }
+}
+
+int fn2ftp(const char *fn)
+{
+ int i, len;
+ const char *feptr;
+ const char *eptr;
+
+ if (!fn)
+ {
+ return efNR;
+ }
+
+ len = std::strlen(fn);
+ if ((len >= 4) && (fn[len - 4] == '.'))
+ {
+ feptr = &(fn[len - 4]);
+ }
+ else
+ {
+ return efNR;
+ }
+
+ for (i = 0; (i < efNR); i++)
+ {
+ if ((eptr = deffile[i].ext) != NULL)
+ {
+ if (gmx_strcasecmp(feptr, eptr) == 0)
+ {
+ break;
+ }
+ }
+ }
+
+ return i;
+}
#include <stdlib.h>
#include <string.h>
+#include <cmath>
+
#include <algorithm>
++#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/domdec.h"
+#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/ewald/ewald.h"
-#include "gromacs/fileio/filenm.h"
-#include "gromacs/gmxlib/gpu_utils/gpu_utils.h"
-#include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/force.h"
-#include "gromacs/legacyheaders/gmx_detect_hardware.h"
-#include "gromacs/legacyheaders/gmx_omp_nthreads.h"
-#include "gromacs/legacyheaders/inputrec.h"
-#include "gromacs/legacyheaders/macros.h"
-#include "gromacs/legacyheaders/md_logging.h"
-#include "gromacs/legacyheaders/md_support.h"
-#include "gromacs/legacyheaders/names.h"
-#include "gromacs/legacyheaders/network.h"
-#include "gromacs/legacyheaders/nonbonded.h"
-#include "gromacs/legacyheaders/ns.h"
-#include "gromacs/legacyheaders/qmmm.h"
-#include "gromacs/legacyheaders/tables.h"
-#include "gromacs/legacyheaders/txtdump.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/fileio/filetypes.h"
+#include "gromacs/gmxlib/md_logging.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/gmxlib/nonbonded/nonbonded.h"
+#include "gromacs/gpu_utils/gpu_utils.h"
+#include "gromacs/hardware/detecthardware.h"
#include "gromacs/listed-forces/manage-threading.h"
+#include "gromacs/listed-forces/pairs.h"
#include "gromacs/math/calculate-ewald-splitting-coefficient.h"
+#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/simd/simd.h"
+#include "gromacs/tables/forcetable.h"
#include "gromacs/topology/mtop_util.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/cstringutil.h"
+ #include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"
const t_commrec *cr,
matrix box,
const char *tabfn,
- const char *tabafn,
const char *tabpfn,
- const char *tabbfn,
+ const t_filenm *tabbfnm,
const char *nbpu_opt,
gmx_bool bNoSolvOpt,
real print_force)
fr->nwall = ir->nwall;
if (ir->nwall && ir->wall_type == ewtTABLE)
{
- make_wall_tables(fp, oenv, ir, tabfn, &mtop->groups, fr);
+ make_wall_tables(fp, ir, tabfn, &mtop->groups, fr);
}
- if (fcd && tabbfn)
+ if (fcd && tabbfnm)
{
- fcd->bondtab = make_bonded_tables(fp,
- F_TABBONDS, F_TABBONDSNC,
- mtop, tabbfn, "b");
- fcd->angletab = make_bonded_tables(fp,
- F_TABANGLES, -1,
- mtop, tabbfn, "a");
- fcd->dihtab = make_bonded_tables(fp,
- F_TABDIHS, -1,
- mtop, tabbfn, "d");
+ // Need to catch std::bad_alloc
+ // TODO Don't need to catch this here, when merging with master branch
+ try
+ {
+ fcd->bondtab = make_bonded_tables(fp,
+ F_TABBONDS, F_TABBONDSNC,
+ mtop, tabbfnm, "b");
+ fcd->angletab = make_bonded_tables(fp,
+ F_TABANGLES, -1,
+ mtop, tabbfnm, "a");
+ fcd->dihtab = make_bonded_tables(fp,
+ F_TABDIHS, -1,
+ mtop, tabbfnm, "d");
+ }
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
}
else
{
--- /dev/null
- * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+/*
+ * 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.
- * \param[in] tabbfn Table potential file for bonded interactions
++ * Copyright (c) 2013,2014,2015,2016, 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.
+ */
+#ifndef GMX_MDLIB_FORCEREC_H
+#define GMX_MDLIB_FORCEREC_H
+
+#include "gromacs/mdlib/force_flags.h"
+#include "gromacs/mdlib/genborn.h"
+#include "gromacs/mdlib/tgroup.h"
+#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/timing/wallcycle.h"
+
+struct t_commrec;
+struct t_fcdata;
++struct t_filenm;
+
+/*! \brief Create a new forcerec structure */
+t_forcerec *mk_forcerec(void);
+
+/*! \brief Print the contents of the forcerec to a file
+ *
+ * \param[in] fplog The log file to print to
+ * \param[in] fr The forcerec structure
+ */
+void pr_forcerec(FILE *fplog, t_forcerec *fr);
+
+/*! \brief Set the number of charge groups and atoms.
+ *
+ * The force calculation needs information on which atoms it
+ * should do work.
+ * \param[inout] fr The forcerec
+ * \param[in] ncg_home Number of charge groups on this processor
+ * \param[in] ncg_force Number of charge groups to compute force on
+ * \param[in] natoms_force Number of atoms to compute force on
+ * \param[in] natoms_force_constr Number of atoms involved in constraints
+ * \param[in] natoms_f_novirsum Number of atoms for which
+ * force is to be compute but no virial
+ */
+void
+forcerec_set_ranges(t_forcerec *fr,
+ int ncg_home, int ncg_force,
+ int natoms_force,
+ int natoms_force_constr, int natoms_f_novirsum);
+
+/*! \brief Initiate table constants
+ *
+ * Initializes the tables in the interaction constant data structure.
+ * \param[in] fp File for debugging output
+ * \param[in] ic Structure holding the table constant
+ * \param[in] rtab The additional distance to add to tables
+ */
+void init_interaction_const_tables(FILE *fp,
+ interaction_const_t *ic,
+ real rtab);
+
+/*! \brief Initialize forcerec structure.
+ *
+ * The Force rec struct must be created with mk_forcerec.
+ * \param[in] fplog File for printing
+ * \param[out] fr The forcerec
+ * \param[in] fcd Force constant data
+ * \param[in] ir Inputrec structure
+ * \param[in] mtop Molecular topology
+ * \param[in] cr Communication structures
+ * \param[in] box Simulation box
+ * \param[in] tabfn Table potential file for non-bonded interactions
+ * \param[in] tabpfn Table potential file for pair interactions
- const char *tabbfn,
++ * \param[in] tabbfnm Table potential files for bonded interactions
+ * \param[in] nbpu_opt Nonbonded Processing Unit (GPU/CPU etc.)
+ * \param[in] bNoSolvOpt Do not use solvent optimization
+ * \param[in] print_force Print forces for atoms with force >= print_force
+ */
+void init_forcerec(FILE *fplog,
+ t_forcerec *fr,
+ t_fcdata *fcd,
+ const t_inputrec *ir,
+ const gmx_mtop_t *mtop,
+ const t_commrec *cr,
+ matrix box,
+ const char *tabfn,
+ const char *tabpfn,
++ const t_filenm *tabbfnm,
+ const char *nbpu_opt,
+ gmx_bool bNoSolvOpt,
+ real print_force);
+
+/*! \brief Divide exclusions over threads
+ *
+ * Set the exclusion load for the local exclusions and possibly threads
+ * \param[out] fr The force record
+ * \param[in] top The topology
+ */
+void forcerec_set_excl_load(t_forcerec *fr,
+ const gmx_localtop_t *top);
+
+/*! \brief Update parameters dependent on box
+ *
+ * Updates parameters in the forcerec that are time dependent
+ * \param[out] fr The force record
+ * \param[in] box The simulation box
+ */
+void update_forcerec(t_forcerec *fr, matrix box);
+
+#endif
--- /dev/null
- * Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
+/*
+ * 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.
- /* Set the mass of completely frozen particles to ALMOST_ZERO iso 0
- * to avoid div by zero in lincs or shake.
- * Note that constraints can still move a partially frozen particle.
++ * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "mdatoms.h"
+
+#include <math.h>
+
+#include "gromacs/math/functions.h"
+#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/qmmm.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/smalloc.h"
+
+#define ALMOST_ZERO 1e-30
+
+t_mdatoms *init_mdatoms(FILE *fp, const gmx_mtop_t *mtop, gmx_bool bFreeEnergy)
+{
+ int a;
+ double tmA, tmB;
+ t_atom *atom;
+ t_mdatoms *md;
+ gmx_mtop_atomloop_all_t aloop;
+
+ snew(md, 1);
+
+ md->nenergrp = mtop->groups.grps[egcENER].nr;
+ md->bVCMgrps = FALSE;
+ tmA = 0.0;
+ tmB = 0.0;
+
+ aloop = gmx_mtop_atomloop_all_init(mtop);
+ while (gmx_mtop_atomloop_all_next(aloop, &a, &atom))
+ {
+ if (ggrpnr(&mtop->groups, egcVCM, a) > 0)
+ {
+ md->bVCMgrps = TRUE;
+ }
+
+ if (bFreeEnergy && PERTURBED(*atom))
+ {
+ md->nPerturbed++;
+ if (atom->mB != atom->m)
+ {
+ md->nMassPerturbed++;
+ }
+ if (atom->qB != atom->q)
+ {
+ md->nChargePerturbed++;
+ }
+ if (atom->typeB != atom->type)
+ {
+ md->nTypePerturbed++;
+ }
+ }
+
+ tmA += atom->m;
+ tmB += atom->mB;
+ }
+
+ md->tmassA = tmA;
+ md->tmassB = tmB;
+
+ if (bFreeEnergy && fp)
+ {
+ fprintf(fp,
+ "There are %d atoms and %d charges for free energy perturbation\n",
+ md->nPerturbed, md->nChargePerturbed);
+ }
+
+ md->bOrires = gmx_mtop_ftype_count(mtop, F_ORIRES);
+
+ return md;
+}
+
+void atoms2md(const gmx_mtop_t *mtop, const t_inputrec *ir,
+ int nindex, const int *index,
+ int homenr,
+ t_mdatoms *md)
+{
+ gmx_bool bLJPME;
+ gmx_mtop_atomlookup_t alook;
+ int i;
+ const t_grpopts *opts;
+ const gmx_groups_t *groups;
+ int nthreads gmx_unused;
+
+ bLJPME = EVDW_PME(ir->vdwtype);
+
+ opts = &ir->opts;
+
+ groups = &mtop->groups;
+
+ /* Index==NULL indicates no DD (unless we have a DD node with no
+ * atoms), so also check for homenr. This should be
+ * signaled properly with an extra parameter or nindex==-1.
+ */
+ if (index == NULL && (homenr > 0))
+ {
+ md->nr = mtop->natoms;
+ }
+ else
+ {
+ md->nr = nindex;
+ }
+
+ if (md->nr > md->nalloc)
+ {
+ md->nalloc = over_alloc_dd(md->nr);
+
+ if (md->nMassPerturbed)
+ {
+ srenew(md->massA, md->nalloc);
+ srenew(md->massB, md->nalloc);
+ }
+ srenew(md->massT, md->nalloc);
+ srenew(md->invmass, md->nalloc);
+ srenew(md->chargeA, md->nalloc);
+ srenew(md->typeA, md->nalloc);
+ if (md->nPerturbed)
+ {
+ srenew(md->chargeB, md->nalloc);
+ srenew(md->typeB, md->nalloc);
+ }
+ if (bLJPME)
+ {
+ srenew(md->sqrt_c6A, md->nalloc);
+ srenew(md->sigmaA, md->nalloc);
+ srenew(md->sigma3A, md->nalloc);
+ if (md->nPerturbed)
+ {
+ srenew(md->sqrt_c6B, md->nalloc);
+ srenew(md->sigmaB, md->nalloc);
+ srenew(md->sigma3B, md->nalloc);
+ }
+ }
+ srenew(md->ptype, md->nalloc);
+ if (opts->ngtc > 1)
+ {
+ srenew(md->cTC, md->nalloc);
+ /* We always copy cTC with domain decomposition */
+ }
+ srenew(md->cENER, md->nalloc);
+ if (opts->ngacc > 1)
+ {
+ srenew(md->cACC, md->nalloc);
+ }
+ if (opts->nFreeze &&
+ (opts->ngfrz > 1 ||
+ opts->nFreeze[0][XX] || opts->nFreeze[0][YY] || opts->nFreeze[0][ZZ]))
+ {
+ srenew(md->cFREEZE, md->nalloc);
+ }
+ if (md->bVCMgrps)
+ {
+ srenew(md->cVCM, md->nalloc);
+ }
+ if (md->bOrires)
+ {
+ srenew(md->cORF, md->nalloc);
+ }
+ if (md->nPerturbed)
+ {
+ srenew(md->bPerturbed, md->nalloc);
+ }
+
+ /* Note that these user t_mdatoms array pointers are NULL
+ * when there is only one group present.
+ * Therefore, when adding code, the user should use something like:
+ * gprnrU1 = (md->cU1==NULL ? 0 : md->cU1[localatindex])
+ */
+ if (mtop->groups.grpnr[egcUser1] != NULL)
+ {
+ srenew(md->cU1, md->nalloc);
+ }
+ if (mtop->groups.grpnr[egcUser2] != NULL)
+ {
+ srenew(md->cU2, md->nalloc);
+ }
+
+ if (ir->bQMMM)
+ {
+ srenew(md->bQM, md->nalloc);
+ }
+ }
+
+ alook = gmx_mtop_atomlookup_init(mtop);
+
+ // cppcheck-suppress unreadVariable
+ nthreads = gmx_omp_nthreads_get(emntDefault);
+#pragma omp parallel for num_threads(nthreads) schedule(static)
+ for (i = 0; i < md->nr; i++)
+ {
+ try
+ {
+ int g, ag;
+ real mA, mB, fac;
+ real c6, c12;
+ t_atom *atom;
+
+ if (index == NULL)
+ {
+ ag = i;
+ }
+ else
+ {
+ ag = index[i];
+ }
+ gmx_mtop_atomnr_to_atom(alook, ag, &atom);
+
+ if (md->cFREEZE)
+ {
+ md->cFREEZE[i] = ggrpnr(groups, egcFREEZE, ag);
+ }
+ if (EI_ENERGY_MINIMIZATION(ir->eI))
+ {
+ /* Displacement is proportional to F, masses used for constraints */
+ mA = 1.0;
+ mB = 1.0;
+ }
+ else if (ir->eI == eiBD)
+ {
+ /* With BD the physical masses are irrelevant.
+ * To keep the code simple we use most of the normal MD code path
+ * for BD. Thus for constraining the masses should be proportional
+ * to the friction coefficient. We set the absolute value such that
+ * m/2<(dx/dt)^2> = m/2*2kT/fric*dt = kT/2 => m=fric*dt/2
+ * Then if we set the (meaningless) velocity to v=dx/dt, we get the
+ * correct kinetic energy and temperature using the usual code path.
+ * Thus with BD v*dt will give the displacement and the reported
+ * temperature can signal bad integration (too large time step).
+ */
+ if (ir->bd_fric > 0)
+ {
+ mA = 0.5*ir->bd_fric*ir->delta_t;
+ mB = 0.5*ir->bd_fric*ir->delta_t;
+ }
+ else
+ {
+ /* The friction coefficient is mass/tau_t */
+ fac = ir->delta_t/opts->tau_t[md->cTC ? groups->grpnr[egcTC][ag] : 0];
+ mA = 0.5*atom->m*fac;
+ mB = 0.5*atom->mB*fac;
+ }
+ }
+ else
+ {
+ mA = atom->m;
+ mB = atom->mB;
+ }
+ if (md->nMassPerturbed)
+ {
+ md->massA[i] = mA;
+ md->massB[i] = mB;
+ }
+ md->massT[i] = mA;
+ if (mA == 0.0)
+ {
+ md->invmass[i] = 0;
+ }
+ else if (md->cFREEZE)
+ {
+ g = md->cFREEZE[i];
+ if (opts->nFreeze[g][XX] && opts->nFreeze[g][YY] && opts->nFreeze[g][ZZ])
+ {
++ /* Set the mass of completely frozen particles to ALMOST_ZERO
++ * iso 0 to avoid div by zero in lincs or shake.
+ */
+ md->invmass[i] = ALMOST_ZERO;
+ }
+ else
+ {
++ /* Note: Partially frozen particles use the normal invmass.
++ * If such particles are constrained, the frozen dimensions
++ * should not be updated with the constrained coordinates.
++ */
+ md->invmass[i] = 1.0/mA;
+ }
+ }
+ else
+ {
+ md->invmass[i] = 1.0/mA;
+ }
+ md->chargeA[i] = atom->q;
+ md->typeA[i] = atom->type;
+ if (bLJPME)
+ {
+ c6 = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c6;
+ c12 = mtop->ffparams.iparams[atom->type*(mtop->ffparams.atnr+1)].lj.c12;
+ md->sqrt_c6A[i] = sqrt(c6);
+ if (c6 == 0.0 || c12 == 0)
+ {
+ md->sigmaA[i] = 1.0;
+ }
+ else
+ {
+ md->sigmaA[i] = gmx::sixthroot(c12/c6);
+ }
+ md->sigma3A[i] = 1/(md->sigmaA[i]*md->sigmaA[i]*md->sigmaA[i]);
+ }
+ if (md->nPerturbed)
+ {
+ md->bPerturbed[i] = PERTURBED(*atom);
+ md->chargeB[i] = atom->qB;
+ md->typeB[i] = atom->typeB;
+ if (bLJPME)
+ {
+ c6 = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c6;
+ c12 = mtop->ffparams.iparams[atom->typeB*(mtop->ffparams.atnr+1)].lj.c12;
+ md->sqrt_c6B[i] = sqrt(c6);
+ if (c6 == 0.0 || c12 == 0)
+ {
+ md->sigmaB[i] = 1.0;
+ }
+ else
+ {
+ md->sigmaB[i] = gmx::sixthroot(c12/c6);
+ }
+ md->sigma3B[i] = 1/(md->sigmaB[i]*md->sigmaB[i]*md->sigmaB[i]);
+ }
+ }
+ md->ptype[i] = atom->ptype;
+ if (md->cTC)
+ {
+ md->cTC[i] = groups->grpnr[egcTC][ag];
+ }
+ md->cENER[i] =
+ (groups->grpnr[egcENER] ? groups->grpnr[egcENER][ag] : 0);
+ if (md->cACC)
+ {
+ md->cACC[i] = groups->grpnr[egcACC][ag];
+ }
+ if (md->cVCM)
+ {
+ md->cVCM[i] = groups->grpnr[egcVCM][ag];
+ }
+ if (md->cORF)
+ {
+ md->cORF[i] = groups->grpnr[egcORFIT][ag];
+ }
+
+ if (md->cU1)
+ {
+ md->cU1[i] = groups->grpnr[egcUser1][ag];
+ }
+ if (md->cU2)
+ {
+ md->cU2[i] = groups->grpnr[egcUser2][ag];
+ }
+
+ if (ir->bQMMM)
+ {
+ if (groups->grpnr[egcQMMM] == 0 ||
+ groups->grpnr[egcQMMM][ag] < groups->grps[egcQMMM].nr-1)
+ {
+ md->bQM[i] = TRUE;
+ }
+ else
+ {
+ md->bQM[i] = FALSE;
+ }
+ }
+ }
+ GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+ }
+
+ gmx_mtop_atomlookup_destroy(alook);
+
+ md->homenr = homenr;
+ md->lambda = 0;
+}
+
+void update_mdatoms(t_mdatoms *md, real lambda)
+{
+ int al, end;
+ real L1 = 1.0-lambda;
+
+ end = md->nr;
+
+ if (md->nMassPerturbed)
+ {
+ for (al = 0; (al < end); al++)
+ {
+ if (md->bPerturbed[al])
+ {
+ md->massT[al] = L1*md->massA[al]+ lambda*md->massB[al];
+ if (md->invmass[al] > 1.1*ALMOST_ZERO)
+ {
+ md->invmass[al] = 1.0/md->massT[al];
+ }
+ }
+ }
+ md->tmass = L1*md->tmassA + lambda*md->tmassB;
+ }
+ else
+ {
+ md->tmass = md->tmassA;
+ }
+ md->lambda = lambda;
+}
/*! Texture reference for Ewald coulomb force table; bound to cu_nbparam_t.coulomb_tab */
texture<float, 1, cudaReadModeElementType> coulomb_tab_texref;
-/* Convenience defines */
-#define NCL_PER_SUPERCL (NBNXN_GPU_NCLUSTER_PER_SUPERCLUSTER)
-#define CL_SIZE (NBNXN_GPU_CLUSTER_SIZE)
-/***** The kernels come here *****/
+/***** The kernel declarations/definitions come here *****/
+ #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_utils.cuh"
-/* Top-level kernel generation: will generate through multiple inclusion the
- * following flavors for all kernels:
+/* Top-level kernel declaration generation: will generate through multiple
+ * inclusion the following flavors for all kernel declarations:
* - force-only output;
* - force and energy output;
* - force-only with pair list pruning;
* shuffle-based reduction, hence CC >= 3.0.
*/
-/* Kernel launch bounds as function of NTHREAD_Z.
+/* Kernel launch bounds for different compute capabilities. The value of NTHREAD_Z
+ * determines the number of threads per block and it is chosen such that
+ * 16 blocks/multiprocessor can be kept in flight.
- * - CC 2.x, 3.0, 3.5, 5.x: NTHREAD_Z=1, (64, 16) bounds
- * - CC 3.7: NTHREAD_Z=2, (128, 16) bounds
+ * - CC 3.0/3.5/5.x, >=6.1: NTHREAD_Z=1, (64, 16) bounds
+ * - CC 3.7, 6.0: NTHREAD_Z=2, (128, 16) bounds
+ *
+ * Note: convenience macros, need to be undef-ed at the end of the file.
*/
- #if GMX_PTX_ARCH == 370
-#if __CUDA_ARCH__ == 370 || __CUDA_ARCH__ == 600
-#define NTHREAD_Z (2)
-#define MIN_BLOCKS_PER_MP (16)
++#if GMX_PTX_ARCH == 370 || GMX_PTX_ARCH == 600
+ #define NTHREAD_Z (2)
+ #define MIN_BLOCKS_PER_MP (16)
#else
-#define NTHREAD_Z (1)
-#define MIN_BLOCKS_PER_MP (16)
-#endif
-#define THREADS_PER_BLOCK (CL_SIZE*CL_SIZE*NTHREAD_Z)
+ #define NTHREAD_Z (1)
+ #define MIN_BLOCKS_PER_MP (16)
- #endif /* GMX_PTX_ARCH == 370 */
++#endif /* GMX_PTX_ARCH == 370 || GMX_PTX_ARCH == 600 */
+#define THREADS_PER_BLOCK (c_clSize*c_clSize*NTHREAD_Z)
-#if __CUDA_ARCH__ >= 350
+#if GMX_PTX_ARCH >= 350
+#if (GMX_PTX_ARCH <= 210) && (NTHREAD_Z > 1)
+ #error NTHREAD_Z > 1 will give incorrect results on CC 2.x
+#endif
+/**@}*/
__launch_bounds__(THREADS_PER_BLOCK, MIN_BLOCKS_PER_MP)
#else
__launch_bounds__(THREADS_PER_BLOCK)
--- /dev/null
- ssd += fabs(2*(f - numf)/(f + numf));
+/*
+ * 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,2015,2016, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "forcetable.h"
+
+#include <cmath>
+
+#include <algorithm>
+
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/nblist.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
+
+/* All the possible (implemented) table functions */
+enum {
+ etabLJ6,
+ etabLJ12,
+ etabLJ6Shift,
+ etabLJ12Shift,
+ etabShift,
+ etabRF,
+ etabRF_ZERO,
+ etabCOUL,
+ etabEwald,
+ etabEwaldSwitch,
+ etabEwaldUser,
+ etabEwaldUserSwitch,
+ etabLJ6Ewald,
+ etabLJ6Switch,
+ etabLJ12Switch,
+ etabCOULSwitch,
+ etabLJ6Encad,
+ etabLJ12Encad,
+ etabCOULEncad,
+ etabEXPMIN,
+ etabUSER,
+ etabNR
+};
+
+/** Evaluates to true if the table type contains user data. */
+#define ETAB_USER(e) ((e) == etabUSER || \
+ (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
+
+typedef struct {
+ const char *name;
+ gmx_bool bCoulomb;
+} t_tab_props;
+
+/* This structure holds name and a flag that tells whether
+ this is a Coulomb type funtion */
+static const t_tab_props tprops[etabNR] = {
+ { "LJ6", FALSE },
+ { "LJ12", FALSE },
+ { "LJ6Shift", FALSE },
+ { "LJ12Shift", FALSE },
+ { "Shift", TRUE },
+ { "RF", TRUE },
+ { "RF-zero", TRUE },
+ { "COUL", TRUE },
+ { "Ewald", TRUE },
+ { "Ewald-Switch", TRUE },
+ { "Ewald-User", TRUE },
+ { "Ewald-User-Switch", TRUE },
+ { "LJ6Ewald", FALSE },
+ { "LJ6Switch", FALSE },
+ { "LJ12Switch", FALSE },
+ { "COULSwitch", TRUE },
+ { "LJ6-Encad shift", FALSE },
+ { "LJ12-Encad shift", FALSE },
+ { "COUL-Encad shift", TRUE },
+ { "EXPMIN", FALSE },
+ { "USER", FALSE },
+};
+
+typedef struct {
+ int nx, nx0;
+ double tabscale;
+ double *x, *v, *f;
+} t_tabledata;
+
+double v_q_ewald_lr(double beta, double r)
+{
+ if (r == 0)
+ {
+ return beta*2/sqrt(M_PI);
+ }
+ else
+ {
+ return std::erf(beta*r)/r;
+ }
+}
+
+double v_lj_ewald_lr(double beta, double r)
+{
+ double br, br2, br4, r6, factor;
+ if (r == 0)
+ {
+ return gmx::power6(beta)/6;
+ }
+ else
+ {
+ br = beta*r;
+ br2 = br*br;
+ br4 = br2*br2;
+ r6 = gmx::power6(r);
+ factor = (1.0 - std::exp(-br2)*(1 + br2 + 0.5*br4))/r6;
+ return factor;
+ }
+}
+
+void table_spline3_fill_ewald_lr(real *table_f,
+ real *table_v,
+ real *table_fdv0,
+ int ntab,
+ double dx,
+ real beta,
+ real_space_grid_contribution_computer v_lr)
+{
+ real tab_max;
+ int i, i_inrange;
+ double dc, dc_new;
+ gmx_bool bOutOfRange;
+ double v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
+ double x_r0;
+
+ /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
+ * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
+ */
+
+ if (ntab < 2)
+ {
+ gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
+ }
+
+ /* We need some margin to be able to divide table values by r
+ * in the kernel and also to do the integration arithmetics
+ * without going out of range. Furthemore, we divide by dx below.
+ */
+ tab_max = GMX_REAL_MAX*0.0001;
+
+ /* This function produces a table with:
+ * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
+ * maximum force error: V'''/(6*4)*dx^2
+ * The rms force error is the max error times 1/sqrt(5)=0.45.
+ */
+
+ bOutOfRange = FALSE;
+ i_inrange = ntab;
+ v_inrange = 0;
+ dc = 0;
+ for (i = ntab-1; i >= 0; i--)
+ {
+ x_r0 = i*dx;
+
+ v_r0 = (*v_lr)(beta, x_r0);
+
+ if (!bOutOfRange)
+ {
+ i_inrange = i;
+ v_inrange = v_r0;
+
+ vi = v_r0;
+ }
+ else
+ {
+ /* Linear continuation for the last point in range */
+ vi = v_inrange - dc*(i - i_inrange)*dx;
+ }
+
+ if (table_v != NULL)
+ {
+ table_v[i] = vi;
+ }
+
+ if (i == 0)
+ {
+ continue;
+ }
+
+ /* Get the potential at table point i-1 */
+ v_r1 = (*v_lr)(beta, (i-1)*dx);
+
+ if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
+ {
+ bOutOfRange = TRUE;
+ }
+
+ if (!bOutOfRange)
+ {
+ /* Calculate the average second derivative times dx over interval i-1 to i.
+ * Using the function values at the end points and in the middle.
+ */
+ a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx);
+ /* Set the derivative of the spline to match the difference in potential
+ * over the interval plus the average effect of the quadratic term.
+ * This is the essential step for minimizing the error in the force.
+ */
+ dc = (v_r0 - v_r1)/dx + 0.5*a2dx;
+ }
+
+ if (i == ntab - 1)
+ {
+ /* Fill the table with the force, minus the derivative of the spline */
+ table_f[i] = -dc;
+ }
+ else
+ {
+ /* tab[i] will contain the average of the splines over the two intervals */
+ table_f[i] += -0.5*dc;
+ }
+
+ if (!bOutOfRange)
+ {
+ /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
+ * matching the potential at the two end points
+ * and the derivative dc at the end point xr.
+ */
+ a0 = v_r0;
+ a1 = dc;
+ a2dx = (a1*dx + v_r1 - a0)*2/dx;
+
+ /* Set dc to the derivative at the next point */
+ dc_new = a1 - a2dx;
+
+ if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
+ {
+ bOutOfRange = TRUE;
+ }
+ else
+ {
+ dc = dc_new;
+ }
+ }
+
+ table_f[(i-1)] = -0.5*dc;
+ }
+ /* Currently the last value only contains half the force: double it */
+ table_f[0] *= 2;
+
+ if (table_v != NULL && table_fdv0 != NULL)
+ {
+ /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
+ * init_ewald_f_table().
+ */
+ for (i = 0; i < ntab-1; i++)
+ {
+ table_fdv0[4*i] = table_f[i];
+ table_fdv0[4*i+1] = table_f[i+1]-table_f[i];
+ table_fdv0[4*i+2] = table_v[i];
+ table_fdv0[4*i+3] = 0.0;
+ }
+ table_fdv0[4*(ntab-1)] = table_f[(ntab-1)];
+ table_fdv0[4*(ntab-1)+1] = -table_f[(ntab-1)];
+ table_fdv0[4*(ntab-1)+2] = table_v[(ntab-1)];
+ table_fdv0[4*(ntab-1)+3] = 0.0;
+ }
+}
+
+/* Returns the spacing for a function using the maximum of
+ * the third derivative, x_scale (unit 1/length)
+ * and function tolerance.
+ */
+static double spline3_table_scale(double third_deriv_max,
+ double x_scale,
+ double func_tol)
+{
+ double deriv_tol;
+ double sc_deriv, sc_func;
+
+ /* Force tolerance: single precision accuracy */
+ deriv_tol = GMX_FLOAT_EPS;
+ sc_deriv = sqrt(third_deriv_max/(6*4*deriv_tol*x_scale))*x_scale;
+
+ /* Don't try to be more accurate on energy than the precision */
+ func_tol = std::max(func_tol, static_cast<double>(GMX_REAL_EPS));
+ sc_func = std::cbrt(third_deriv_max/(6*12*sqrt(3)*func_tol))*x_scale;
+
+ return std::max(sc_deriv, sc_func);
+}
+
+/* The scale (1/spacing) for third order spline interpolation
+ * of the Ewald mesh contribution which needs to be subtracted
+ * from the non-bonded interactions.
+ * Since there is currently only one spacing for Coulomb and LJ,
+ * the finest spacing is used if both Ewald types are present.
+ *
+ * Note that we could also implement completely separate tables
+ * for Coulomb and LJ Ewald, each with their own spacing.
+ * The current setup with the same spacing can provide slightly
+ * faster kernels with both Coulomb and LJ Ewald, especially
+ * when interleaving both tables (currently not implemented).
+ */
+real ewald_spline3_table_scale(const interaction_const_t *ic)
+{
+ real sc;
+
+ sc = 0;
+
+ if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype))
+ {
+ double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
+ double etol;
+ real sc_q;
+
+ /* Energy tolerance: 0.1 times the cut-off jump */
+ etol = 0.1*std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+
+ sc_q = spline3_table_scale(erf_x_d3, ic->ewaldcoeff_q, etol);
+
+ if (debug)
+ {
+ fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f 1/nm\n", 1/sc_q);
+ }
+
+ sc = std::max(sc, sc_q);
+ }
+
+ if (EVDW_PME(ic->vdwtype))
+ {
+ double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
+ double xrc2, etol;
+ real sc_lj;
+
+ /* Energy tolerance: 0.1 times the cut-off jump */
+ xrc2 = gmx::square(ic->ewaldcoeff_lj*ic->rvdw);
+ etol = 0.1*std::exp(-xrc2)*(1 + xrc2 + xrc2*xrc2/2.0);
+
+ sc_lj = spline3_table_scale(func_d3, ic->ewaldcoeff_lj, etol);
+
+ if (debug)
+ {
+ fprintf(debug, "Ewald LJ quadratic spline table spacing: %f 1/nm\n", 1/sc_lj);
+ }
+
+ sc = std::max(sc, sc_lj);
+ }
+
+ return sc;
+}
+
+/* Calculate the potential and force for an r value
+ * in exactly the same way it is done in the inner loop.
+ * VFtab is a pointer to the table data, offset is
+ * the point where we should begin and stride is
+ * 4 if we have a buckingham table, 3 otherwise.
+ * If you want to evaluate table no N, set offset to 4*N.
+ *
+ * We use normal precision here, since that is what we
+ * will use in the inner loops.
+ */
+static void evaluate_table(real VFtab[], int offset, int stride,
+ real tabscale, real r, real *y, real *yp)
+{
+ int n;
+ real rt, eps, eps2;
+ real Y, F, Geps, Heps2, Fp;
+
+ rt = r*tabscale;
+ n = (int)rt;
+ eps = rt - n;
+ eps2 = eps*eps;
+ n = offset+stride*n;
+ Y = VFtab[n];
+ F = VFtab[n+1];
+ Geps = eps*VFtab[n+2];
+ Heps2 = eps2*VFtab[n+3];
+ Fp = F+Geps+Heps2;
+ *y = Y+eps*Fp;
+ *yp = (Fp+Geps+2.0*Heps2)*tabscale;
+}
+
+static void copy2table(int n, int offset, int stride,
+ double x[], double Vtab[], double Ftab[], real scalefactor,
+ real dest[])
+{
+/* Use double prec. for the intermediary variables
+ * and temporary x/vtab/vtab2 data to avoid unnecessary
+ * loss of precision.
+ */
+ int i, nn0;
+ double F, G, H, h;
+
+ h = 0;
+ for (i = 0; (i < n); i++)
+ {
+ if (i < n-1)
+ {
+ h = x[i+1] - x[i];
+ F = -Ftab[i]*h;
+ G = 3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
+ H = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] + Ftab[i])*h;
+ }
+ else
+ {
+ /* Fill the last entry with a linear potential,
+ * this is mainly for rounding issues with angle and dihedral potentials.
+ */
+ F = -Ftab[i]*h;
+ G = 0;
+ H = 0;
+ }
+ nn0 = offset + i*stride;
+ dest[nn0] = scalefactor*Vtab[i];
+ dest[nn0+1] = scalefactor*F;
+ dest[nn0+2] = scalefactor*G;
+ dest[nn0+3] = scalefactor*H;
+ }
+}
+
+static void init_table(int n, int nx0,
+ double tabscale, t_tabledata *td, gmx_bool bAlloc)
+{
+ td->nx = n;
+ td->nx0 = nx0;
+ td->tabscale = tabscale;
+ if (bAlloc)
+ {
+ snew(td->x, td->nx);
+ snew(td->v, td->nx);
+ snew(td->f, td->nx);
+ }
+}
+
+static void spline_forces(int nx, double h, double v[], gmx_bool bS3, gmx_bool bE3,
+ double f[])
+{
+ int start, end, i;
+ double v3, b_s, b_e, b;
+ double beta, *gamma;
+
+ /* Formulas can be found in:
+ * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
+ */
+
+ if (nx < 4 && (bS3 || bE3))
+ {
+ gmx_fatal(FARGS, "Can not generate splines with third derivative boundary conditions with less than 4 (%d) points", nx);
+ }
+
+ /* To make life easy we initially set the spacing to 1
+ * and correct for this at the end.
+ */
+ if (bS3)
+ {
+ /* Fit V''' at the start */
+ v3 = v[3] - 3*v[2] + 3*v[1] - v[0];
+ if (debug)
+ {
+ fprintf(debug, "The left third derivative is %g\n", v3/(h*h*h));
+ }
+ b_s = 2*(v[1] - v[0]) + v3/6;
+ start = 0;
+
+ if (FALSE)
+ {
+ /* Fit V'' at the start */
+ real v2;
+
+ v2 = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
+ /* v2 = v[2] - 2*v[1] + v[0]; */
+ if (debug)
+ {
+ fprintf(debug, "The left second derivative is %g\n", v2/(h*h));
+ }
+ b_s = 3*(v[1] - v[0]) - v2/2;
+ start = 0;
+ }
+ }
+ else
+ {
+ b_s = 3*(v[2] - v[0]) + f[0]*h;
+ start = 1;
+ }
+ if (bE3)
+ {
+ /* Fit V''' at the end */
+ v3 = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
+ if (debug)
+ {
+ fprintf(debug, "The right third derivative is %g\n", v3/(h*h*h));
+ }
+ b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
+ end = nx;
+ }
+ else
+ {
+ /* V'=0 at the end */
+ b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
+ end = nx - 1;
+ }
+
+ snew(gamma, nx);
+ beta = (bS3 ? 1 : 4);
+
+ /* For V'' fitting */
+ /* beta = (bS3 ? 2 : 4); */
+
+ f[start] = b_s/beta;
+ for (i = start+1; i < end; i++)
+ {
+ gamma[i] = 1/beta;
+ beta = 4 - gamma[i];
+ b = 3*(v[i+1] - v[i-1]);
+ f[i] = (b - f[i-1])/beta;
+ }
+ gamma[end-1] = 1/beta;
+ beta = (bE3 ? 1 : 4) - gamma[end-1];
+ f[end-1] = (b_e - f[end-2])/beta;
+
+ for (i = end-2; i >= start; i--)
+ {
+ f[i] -= gamma[i+1]*f[i+1];
+ }
+ sfree(gamma);
+
+ /* Correct for the minus sign and the spacing */
+ for (i = start; i < end; i++)
+ {
+ f[i] = -f[i]/h;
+ }
+}
+
+static void set_forces(FILE *fp, int angle,
+ int nx, double h, double v[], double f[],
+ int table)
+{
+ int start, end;
+
+ if (angle == 2)
+ {
+ gmx_fatal(FARGS,
+ "Force generation for dihedral tables is not (yet) implemented");
+ }
+
+ start = 0;
+ while (v[start] == 0)
+ {
+ start++;
+ }
+
+ end = nx;
+ while (v[end-1] == 0)
+ {
+ end--;
+ }
+ if (end > nx - 2)
+ {
+ end = nx;
+ }
+ else
+ {
+ end++;
+ }
+
+ if (fp)
+ {
+ fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
+ table+1, start*h, end == nx ? "V'''" : "V'=0", (end-1)*h);
+ }
+ spline_forces(end-start, h, v+start, TRUE, end == nx, f+start);
+}
+
+static void read_tables(FILE *fp, const char *fn,
+ int ntab, int angle, t_tabledata td[])
+{
+ char *libfn;
+ char buf[STRLEN];
+ double **yy = NULL, start, end, dx0, dx1, ssd, vm, vp, f, numf;
+ int k, i, nx, nx0 = 0, ny, nny, ns;
+ gmx_bool bAllZero, bZeroV, bZeroF;
+ double tabscale;
+
+ nny = 2*ntab+1;
+ libfn = gmxlibfn(fn);
+ nx = read_xvg(libfn, &yy, &ny);
+ if (ny != nny)
+ {
+ gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d",
+ libfn, ny, nny);
+ }
+ if (angle == 0)
+ {
+ if (yy[0][0] != 0.0)
+ {
+ gmx_fatal(FARGS,
+ "The first distance in file %s is %f nm instead of %f nm",
+ libfn, yy[0][0], 0.0);
+ }
+ }
+ else
+ {
+ if (angle == 1)
+ {
+ start = 0.0;
+ }
+ else
+ {
+ start = -180.0;
+ }
+ end = 180.0;
+ if (yy[0][0] != start || yy[0][nx-1] != end)
+ {
+ gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
+ libfn, start, end, yy[0][0], yy[0][nx-1]);
+ }
+ }
+
+ tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
+
+ if (fp)
+ {
+ fprintf(fp, "Read user tables from %s with %d data points.\n", libfn, nx);
+ if (angle == 0)
+ {
+ fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
+ }
+ }
+
+ bAllZero = TRUE;
+ for (k = 0; k < ntab; k++)
+ {
+ bZeroV = TRUE;
+ bZeroF = TRUE;
+ for (i = 0; (i < nx); i++)
+ {
+ if (i >= 2)
+ {
+ dx0 = yy[0][i-1] - yy[0][i-2];
+ dx1 = yy[0][i] - yy[0][i-1];
+ /* Check for 1% deviation in spacing */
+ if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1)))
+ {
+ gmx_fatal(FARGS, "In table file '%s' the x values are not equally spaced: %f %f %f", fn, yy[0][i-2], yy[0][i-1], yy[0][i]);
+ }
+ }
+ if (yy[1+k*2][i] != 0)
+ {
+ bZeroV = FALSE;
+ if (bAllZero)
+ {
+ bAllZero = FALSE;
+ nx0 = i;
+ }
+ if (yy[1+k*2][i] > 0.01*GMX_REAL_MAX ||
+ yy[1+k*2][i] < -0.01*GMX_REAL_MAX)
+ {
+ gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
+ yy[1+k*2][i], fn);
+ }
+ }
+ if (yy[1+k*2+1][i] != 0)
+ {
+ bZeroF = FALSE;
+ if (bAllZero)
+ {
+ bAllZero = FALSE;
+ nx0 = i;
+ }
+ if (yy[1+k*2+1][i] > 0.01*GMX_REAL_MAX ||
+ yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX)
+ {
+ gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
+ yy[1+k*2+1][i], fn);
+ }
+ }
+ }
+
+ if (!bZeroV && bZeroF)
+ {
+ set_forces(fp, angle, nx, 1/tabscale, yy[1+k*2], yy[1+k*2+1], k);
+ }
+ else
+ {
+ /* Check if the second column is close to minus the numerical
+ * derivative of the first column.
+ */
+ ssd = 0;
+ ns = 0;
+ for (i = 1; (i < nx-1); i++)
+ {
+ vm = yy[1+2*k][i-1];
+ vp = yy[1+2*k][i+1];
+ f = yy[1+2*k+1][i];
+ if (vm != 0 && vp != 0 && f != 0)
+ {
+ /* Take the centered difference */
+ numf = -(vp - vm)*0.5*tabscale;
- sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %d%% from minus the numerical derivative of the potential\n", ns, k, libfn, (int)(100*ssd+0.5));
++ if (f + numf != 0)
++ {
++ ssd += fabs(2*(f - numf)/(f + numf));
++ }
+ ns++;
+ }
+ }
+ if (ns > 0)
+ {
+ ssd /= ns;
- bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle)
++ sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %lld%% from minus the numerical derivative of the potential\n", ns, k, libfn, (long long int)(100*ssd+0.5));
+ if (debug)
+ {
+ fprintf(debug, "%s", buf);
+ }
+ if (ssd > 0.2)
+ {
+ if (fp)
+ {
+ fprintf(fp, "\nWARNING: %s\n", buf);
+ }
+ fprintf(stderr, "\nWARNING: %s\n", buf);
+ }
+ }
+ }
+ }
+ if (bAllZero && fp)
+ {
+ fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn);
+ }
+
+ for (k = 0; (k < ntab); k++)
+ {
+ init_table(nx, nx0, tabscale, &(td[k]), TRUE);
+ for (i = 0; (i < nx); i++)
+ {
+ td[k].x[i] = yy[0][i];
+ td[k].v[i] = yy[2*k+1][i];
+ td[k].f[i] = yy[2*k+2][i];
+ }
+ }
+ for (i = 0; (i < ny); i++)
+ {
+ sfree(yy[i]);
+ }
+ sfree(yy);
+ sfree(libfn);
+}
+
+static void done_tabledata(t_tabledata *td)
+{
+ if (!td)
+ {
+ return;
+ }
+
+ sfree(td->x);
+ sfree(td->v);
+ sfree(td->f);
+}
+
+static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
+ gmx_bool b14only)
+{
+ /* Fill the table according to the formulas in the manual.
+ * In principle, we only need the potential and the second
+ * derivative, but then we would have to do lots of calculations
+ * in the inner loop. By precalculating some terms (see manual)
+ * we get better eventual performance, despite a larger table.
+ *
+ * Since some of these higher-order terms are very small,
+ * we always use double precision to calculate them here, in order
+ * to avoid unnecessary loss of precision.
+ */
+ int i;
+ double reppow, p;
+ double r1, rc, r12, r13;
+ double r, r2, r6, rc2, rc6, rc12;
+ double expr, Vtab, Ftab;
+ /* Parameters for David's function */
+ double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
+ /* Parameters for the switching function */
+ double ksw, swi, swi1;
+ /* Temporary parameters */
+ gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
+ double ewc = fr->ewaldcoeff_q;
+ double ewclj = fr->ewaldcoeff_lj;
+ double Vcut = 0;
+
+ if (b14only)
+ {
+ bPotentialSwitch = FALSE;
+ bForceSwitch = FALSE;
+ bPotentialShift = FALSE;
+ }
+ else
+ {
+ bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
+ (tp == etabCOULSwitch) ||
+ (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
+ bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
+ (tp == etabShift) ||
+ (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
+ bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
+ (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
+ }
+
+ reppow = fr->reppow;
+
+ if (tprops[tp].bCoulomb)
+ {
+ r1 = fr->rcoulomb_switch;
+ rc = fr->rcoulomb;
+ }
+ else
+ {
+ r1 = fr->rvdw_switch;
+ rc = fr->rvdw;
+ }
+ if (bPotentialSwitch)
+ {
+ ksw = 1.0/(gmx::power5(rc-r1));
+ }
+ else
+ {
+ ksw = 0.0;
+ }
+ if (bForceSwitch)
+ {
+ if (tp == etabShift)
+ {
+ p = 1;
+ }
+ else if (tp == etabLJ6Shift)
+ {
+ p = 6;
+ }
+ else
+ {
+ p = reppow;
+ }
+
+ A = p * ((p+1)*r1-(p+4)*rc)/(std::pow(rc, p+2)*gmx::square(rc-r1));
+ B = -p * ((p+1)*r1-(p+3)*rc)/(std::pow(rc, p+2)*gmx::power3(rc-r1));
+ C = 1.0/std::pow(rc, p)-A/3.0*gmx::power3(rc-r1)-B/4.0*gmx::power4(rc-r1);
+ if (tp == etabLJ6Shift)
+ {
+ A = -A;
+ B = -B;
+ C = -C;
+ }
+ A_3 = A/3.0;
+ B_4 = B/4.0;
+ }
+ if (debug)
+ {
+ fprintf(debug, "Setting up tables\n"); fflush(debug);
+ }
+
+ if (bPotentialShift)
+ {
+ rc2 = rc*rc;
+ rc6 = 1.0/(rc2*rc2*rc2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ rc12 = rc6*rc6;
+ }
+ else
+ {
+ rc12 = std::pow(rc, -reppow);
+ }
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vcut = -rc6;
+ break;
+ case etabLJ6Ewald:
+ Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + gmx::power4(ewclj)*rc2*rc2/2);
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vcut = rc12;
+ break;
+ case etabCOUL:
+ Vcut = 1.0/rc;
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vcut = std::erfc(ewc*rc)/rc;
+ break;
+ case etabEwaldUser:
+ /* Only calculate minus the reciprocal space contribution */
+ Vcut = -std::erf(ewc*rc)/rc;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ /* No need for preventing the usage of modifiers with RF */
+ Vcut = 0.0;
+ break;
+ case etabEXPMIN:
+ Vcut = exp(-rc);
+ break;
+ default:
+ gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
+ tprops[tp].name, __FILE__, __LINE__);
+ }
+ }
+
+ for (i = 0; (i < td->nx); i++)
+ {
+ td->x[i] = i/td->tabscale;
+ }
+ for (i = td->nx0; (i < td->nx); i++)
+ {
+ r = td->x[i];
+ r2 = r*r;
+ r6 = 1.0/(r2*r2*r2);
+ if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
+ {
+ r12 = r6*r6;
+ }
+ else
+ {
+ r12 = std::pow(r, -reppow);
+ }
+ Vtab = 0.0;
+ Ftab = 0.0;
+ if (bPotentialSwitch)
+ {
+ /* swi is function, swi1 1st derivative and swi2 2nd derivative */
+ /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
+ * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
+ * r1 and rc.
+ * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
+ */
+ if (r <= r1)
+ {
+ swi = 1.0;
+ swi1 = 0.0;
+ }
+ else if (r >= rc)
+ {
+ swi = 0.0;
+ swi1 = 0.0;
+ }
+ else
+ {
+ swi = 1 - 10*gmx::power3(r-r1)*ksw*gmx::square(rc-r1)
+ + 15*gmx::power4(r-r1)*ksw*(rc-r1) - 6*gmx::power5(r-r1)*ksw;
+ swi1 = -30*gmx::square(r-r1)*ksw*gmx::square(rc-r1)
+ + 60*gmx::power3(r-r1)*ksw*(rc-r1) - 30*gmx::power4(r-r1)*ksw;
+ }
+ }
+ else /* not really needed, but avoids compiler warnings... */
+ {
+ swi = 1.0;
+ swi1 = 0.0;
+ }
+
+ rc6 = rc*rc*rc;
+ rc6 = 1.0/(rc6*rc6);
+
+ switch (tp)
+ {
+ case etabLJ6:
+ /* Dispersion */
+ Vtab = -r6;
+ Ftab = 6.0*Vtab/r;
+ break;
+ case etabLJ6Switch:
+ case etabLJ6Shift:
+ /* Dispersion */
+ if (r < rc)
+ {
+ Vtab = -r6;
+ Ftab = 6.0*Vtab/r;
+ break;
+ }
+ break;
+ case etabLJ12:
+ /* Repulsion */
+ Vtab = r12;
+ Ftab = reppow*Vtab/r;
+ break;
+ case etabLJ12Switch:
+ case etabLJ12Shift:
+ /* Repulsion */
+ if (r < rc)
+ {
+ Vtab = r12;
+ Ftab = reppow*Vtab/r;
+ }
+ break;
+ case etabLJ6Encad:
+ if (r < rc)
+ {
+ Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
+ Ftab = -(6.0*r6/r-6.0*rc6/rc);
+ }
+ else /* r>rc */
+ {
+ Vtab = 0;
+ Ftab = 0;
+ }
+ break;
+ case etabLJ12Encad:
+ if (r < rc)
+ {
+ Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
+ Ftab = -(6.0*r6/r-6.0*rc6/rc);
+ }
+ else /* r>rc */
+ {
+ Vtab = 0;
+ Ftab = 0;
+ }
+ break;
+ case etabCOUL:
+ Vtab = 1.0/r;
+ Ftab = 1.0/r2;
+ break;
+ case etabCOULSwitch:
+ case etabShift:
+ if (r < rc)
+ {
+ Vtab = 1.0/r;
+ Ftab = 1.0/r2;
+ }
+ break;
+ case etabEwald:
+ case etabEwaldSwitch:
+ Vtab = std::erfc(ewc*r)/r;
+ Ftab = std::erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
+ break;
+ case etabEwaldUser:
+ case etabEwaldUserSwitch:
+ /* Only calculate the negative of the reciprocal space contribution */
+ Vtab = -std::erf(ewc*r)/r;
+ Ftab = -std::erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
+ break;
+ case etabLJ6Ewald:
+ Vtab = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + gmx::power4(ewclj)*r2*r2/2);
+ Ftab = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*gmx::power5(ewclj)*ewclj*r2*r2*r;
+ break;
+ case etabRF:
+ case etabRF_ZERO:
+ Vtab = 1.0/r + fr->k_rf*r2 - fr->c_rf;
+ Ftab = 1.0/r2 - 2*fr->k_rf*r;
+ if (tp == etabRF_ZERO && r >= rc)
+ {
+ Vtab = 0;
+ Ftab = 0;
+ }
+ break;
+ case etabEXPMIN:
+ expr = exp(-r);
+ Vtab = expr;
+ Ftab = expr;
+ break;
+ case etabCOULEncad:
+ if (r < rc)
+ {
+ Vtab = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
+ Ftab = 1.0/r2-1.0/(rc*rc);
+ }
+ else /* r>rc */
+ {
+ Vtab = 0;
+ Ftab = 0;
+ }
+ break;
+ default:
+ gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
+ tp, __FILE__, __LINE__);
+ }
+ if (bForceSwitch)
+ {
+ /* Normal coulomb with cut-off correction for potential */
+ if (r < rc)
+ {
+ Vtab -= C;
+ /* If in Shifting range add something to it */
+ if (r > r1)
+ {
+ r12 = (r-r1)*(r-r1);
+ r13 = (r-r1)*r12;
+ Vtab += -A_3*r13 - B_4*r12*r12;
+ Ftab += A*r12 + B*r13;
+ }
+ }
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ }
+ if (bPotentialShift)
+ {
+ if (r < rc)
+ {
+ Vtab -= Vcut;
+ }
+ else
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ }
+
+ if (ETAB_USER(tp))
+ {
+ Vtab += td->v[i];
+ Ftab += td->f[i];
+ }
+
+ if (bPotentialSwitch)
+ {
+ if (r >= rc)
+ {
+ /* Make sure interactions are zero outside cutoff with modifiers */
+ Vtab = 0;
+ Ftab = 0;
+ }
+ else if (r > r1)
+ {
+ Ftab = Ftab*swi - Vtab*swi1;
+ Vtab = Vtab*swi;
+ }
+ }
+ /* Convert to single precision when we store to mem */
+ td->v[i] = Vtab;
+ td->f[i] = Ftab;
+ }
+
+ /* Continue the table linearly from nx0 to 0.
+ * These values are only required for energy minimization with overlap or TPI.
+ */
+ for (i = td->nx0-1; i >= 0; i--)
+ {
+ td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
+ td->f[i] = td->f[i+1];
+ }
+}
+
+static void set_table_type(int tabsel[], const t_forcerec *fr, gmx_bool b14only)
+{
+ int eltype, vdwtype;
+
+ /* Set the different table indices.
+ * Coulomb first.
+ */
+
+
+ if (b14only)
+ {
+ switch (fr->eeltype)
+ {
+ case eelUSER:
+ case eelPMEUSER:
+ case eelPMEUSERSWITCH:
+ eltype = eelUSER;
+ break;
+ default:
+ eltype = eelCUT;
+ }
+ }
+ else
+ {
+ eltype = fr->eeltype;
+ }
+
+ switch (eltype)
+ {
+ case eelCUT:
+ tabsel[etiCOUL] = etabCOUL;
+ break;
+ case eelPOISSON:
+ tabsel[etiCOUL] = etabShift;
+ break;
+ case eelSHIFT:
+ if (fr->rcoulomb > fr->rcoulomb_switch)
+ {
+ tabsel[etiCOUL] = etabShift;
+ }
+ else
+ {
+ tabsel[etiCOUL] = etabCOUL;
+ }
+ break;
+ case eelEWALD:
+ case eelPME:
+ case eelP3M_AD:
+ tabsel[etiCOUL] = etabEwald;
+ break;
+ case eelPMESWITCH:
+ tabsel[etiCOUL] = etabEwaldSwitch;
+ break;
+ case eelPMEUSER:
+ tabsel[etiCOUL] = etabEwaldUser;
+ break;
+ case eelPMEUSERSWITCH:
+ tabsel[etiCOUL] = etabEwaldUserSwitch;
+ break;
+ case eelRF:
+ case eelGRF:
+ case eelRF_ZERO:
+ tabsel[etiCOUL] = etabRF_ZERO;
+ break;
+ case eelSWITCH:
+ tabsel[etiCOUL] = etabCOULSwitch;
+ break;
+ case eelUSER:
+ tabsel[etiCOUL] = etabUSER;
+ break;
+ case eelENCADSHIFT:
+ tabsel[etiCOUL] = etabCOULEncad;
+ break;
+ default:
+ gmx_fatal(FARGS, "Invalid eeltype %d", eltype);
+ }
+
+ /* Van der Waals time */
+ if (fr->bBHAM && !b14only)
+ {
+ tabsel[etiLJ6] = etabLJ6;
+ tabsel[etiLJ12] = etabEXPMIN;
+ }
+ else
+ {
+ if (b14only && fr->vdwtype != evdwUSER)
+ {
+ vdwtype = evdwCUT;
+ }
+ else
+ {
+ vdwtype = fr->vdwtype;
+ }
+
+ switch (vdwtype)
+ {
+ case evdwSWITCH:
+ tabsel[etiLJ6] = etabLJ6Switch;
+ tabsel[etiLJ12] = etabLJ12Switch;
+ break;
+ case evdwSHIFT:
+ tabsel[etiLJ6] = etabLJ6Shift;
+ tabsel[etiLJ12] = etabLJ12Shift;
+ break;
+ case evdwUSER:
+ tabsel[etiLJ6] = etabUSER;
+ tabsel[etiLJ12] = etabUSER;
+ break;
+ case evdwCUT:
+ tabsel[etiLJ6] = etabLJ6;
+ tabsel[etiLJ12] = etabLJ12;
+ break;
+ case evdwENCADSHIFT:
+ tabsel[etiLJ6] = etabLJ6Encad;
+ tabsel[etiLJ12] = etabLJ12Encad;
+ break;
+ case evdwPME:
+ tabsel[etiLJ6] = etabLJ6Ewald;
+ tabsel[etiLJ12] = etabLJ12;
+ break;
+ default:
+ gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype,
+ __FILE__, __LINE__);
+ }
+
+ if (!b14only && fr->vdw_modifier != eintmodNONE)
+ {
+ if (fr->vdw_modifier != eintmodPOTSHIFT &&
+ fr->vdwtype != evdwCUT)
+ {
+ gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
+ }
+
+ /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
+ * to the original interaction forms when we fill the table, so we only check cutoffs here.
+ */
+ if (fr->vdwtype == evdwCUT)
+ {
+ switch (fr->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ case eintmodEXACTCUTOFF:
+ /* No modification */
+ break;
+ case eintmodPOTSWITCH:
+ tabsel[etiLJ6] = etabLJ6Switch;
+ tabsel[etiLJ12] = etabLJ12Switch;
+ break;
+ case eintmodFORCESWITCH:
+ tabsel[etiLJ6] = etabLJ6Shift;
+ tabsel[etiLJ12] = etabLJ12Shift;
+ break;
+ default:
+ gmx_incons("Unsupported vdw_modifier");
+ }
+ }
+ }
+ }
+}
+
+t_forcetable *make_tables(FILE *out,
+ const t_forcerec *fr,
+ const char *fn,
+ real rtab, int flags)
+{
+ t_tabledata *td;
+ gmx_bool b14only, useUserTable;
+ int nx0, tabsel[etiNR];
+ real scalefactor;
+
+ t_forcetable *table;
+
+ snew(table, 1);
+ b14only = (flags & GMX_MAKETABLES_14ONLY);
+
+ if (flags & GMX_MAKETABLES_FORCEUSER)
+ {
+ tabsel[etiCOUL] = etabUSER;
+ tabsel[etiLJ6] = etabUSER;
+ tabsel[etiLJ12] = etabUSER;
+ }
+ else
+ {
+ set_table_type(tabsel, fr, b14only);
+ }
+ snew(td, etiNR);
+ table->r = rtab;
+ table->scale = 0;
+ table->n = 0;
+
+ table->interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
+ table->format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
+ table->formatsize = 4;
+ table->ninteractions = etiNR;
+ table->stride = table->formatsize*table->ninteractions;
+
+ /* Check whether we have to read or generate */
+ useUserTable = FALSE;
+ for (unsigned int i = 0; (i < etiNR); i++)
+ {
+ if (ETAB_USER(tabsel[i]))
+ {
+ useUserTable = TRUE;
+ }
+ }
+ if (useUserTable)
+ {
+ read_tables(out, fn, etiNR, 0, td);
+ if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
+ {
+ table->n = td[0].nx;
+ }
+ else
+ {
+ if (td[0].x[td[0].nx-1] < rtab)
+ {
+ gmx_fatal(FARGS, "Tables in file %s not long enough for cut-off:\n"
+ "\tshould be at least %f nm\n", fn, rtab);
+ }
+ table->n = (int)(rtab*td[0].tabscale + 0.5);
+ }
+ table->scale = td[0].tabscale;
+ nx0 = td[0].nx0;
+ }
+ else
+ {
+ // No tables are read
+#if GMX_DOUBLE
+ table->scale = 2000.0;
+#else
+ table->scale = 500.0;
+#endif
+ table->n = static_cast<int>(rtab*table->scale);
+ nx0 = 10;
+ }
+
+ /* Each table type (e.g. coul,lj6,lj12) requires four
+ * numbers per table->n+1 data points. For performance reasons we want
+ * the table data to be aligned to a 32-byte boundary.
+ */
+ snew_aligned(table->data, table->stride*(table->n+1)*sizeof(real), 32);
+
+ for (int k = 0; (k < etiNR); k++)
+ {
+ /* Now fill data for tables that should not be read (perhaps
+ overwriting data that was read but should not be used) */
+ if (!ETAB_USER(tabsel[k]))
+ {
+ real scale = table->scale;
+ if (fr->bBHAM && (fr->bham_b_max != 0) && tabsel[k] == etabEXPMIN)
+ {
+ scale /= fr->bham_b_max;
+ }
+ init_table(table->n, nx0, scale, &(td[k]), !useUserTable);
+
+ fill_table(&(td[k]), tabsel[k], fr, b14only);
+ if (out)
+ {
+ fprintf(out, "Generated table with %d data points for %s%s.\n"
+ "Tabscale = %g points/nm\n",
+ td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name,
+ td[k].tabscale);
+ }
+ }
+
+ /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
+ * by including the derivative constants (6.0 or 12.0) in the parameters, since
+ * we no longer calculate force in most steps. This means the c6/c12 parameters
+ * have been scaled up, so we need to scale down the table interactions too.
+ * It comes here since we need to scale user tables too.
+ */
+ if (k == etiLJ6)
+ {
+ scalefactor = 1.0/6.0;
+ }
+ else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
+ {
+ scalefactor = 1.0/12.0;
+ }
+ else
+ {
+ scalefactor = 1.0;
+ }
+
+ copy2table(table->n, k*table->formatsize, table->stride, td[k].x, td[k].v, td[k].f, scalefactor, table->data);
+
+ done_tabledata(&(td[k]));
+ }
+ sfree(td);
+
+ return table;
+}
+
+t_forcetable *make_gb_table(const t_forcerec *fr)
+{
+ t_tabledata *td;
+ int nx0;
+ double r, r2, Vtab, Ftab, expterm;
+
+ t_forcetable *table;
+
+ /* Set the table dimensions for GB, not really necessary to
+ * use etiNR (since we only have one table, but ...)
+ */
+ snew(table, 1);
+ snew(td, 1);
+ table->interaction = GMX_TABLE_INTERACTION_ELEC;
+ table->format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
+ table->r = fr->gbtabr;
+ table->scale = fr->gbtabscale;
+ table->n = static_cast<int>(table->scale*table->r);
+ table->formatsize = 4;
+ table->ninteractions = 1;
+ table->stride = table->formatsize*table->ninteractions;
+ nx0 = 0;
+
+ /* Each table type (e.g. coul,lj6,lj12) requires four numbers per
+ * datapoint. For performance reasons we want the table data to be
+ * aligned on a 32-byte boundary. This new pointer must not be
+ * used in a free() call, but thankfully we're sloppy enough not
+ * to do this :-)
+ */
+
+ snew_aligned(table->data, table->stride*table->n, 32);
+
+ init_table(table->n, nx0, table->scale, &(td[0]), TRUE);
+
+ /* Local implementation so we don't have to use the etabGB
+ * enum above, which will cause problems later when
+ * making the other tables (right now even though we are using
+ * GB, the normal Coulomb tables will be created, but this
+ * will cause a problem since fr->eeltype==etabGB which will not
+ * be defined in fill_table and set_table_type
+ */
+
+ for (int i = nx0; i < table->n; i++)
+ {
+ r = td->x[i];
+ r2 = r*r;
+ expterm = exp(-0.25*r2);
+
+ Vtab = 1/sqrt(r2+expterm);
+ Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
+
+ /* Convert to single precision when we store to mem */
+ td->x[i] = i/table->scale;
+ td->v[i] = Vtab;
+ td->f[i] = Ftab;
+
+ }
+
+ copy2table(table->n, 0, table->stride, td[0].x, td[0].v, td[0].f, 1.0, table->data);
+
+ done_tabledata(&(td[0]));
+ sfree(td);
+
+ return table;
+
+
+}
+
++bondedtable_t make_bonded_table(FILE *fplog, const char *fn, int angle)
+{
+ t_tabledata td;
+ int i;
+ bondedtable_t tab;
+ int stride = 4;
+
+ read_tables(fplog, fn, 1, angle, &td);
+ if (angle > 0)
+ {
+ /* Convert the table from degrees to radians */
+ for (i = 0; i < td.nx; i++)
+ {
+ td.x[i] *= DEG2RAD;
+ td.f[i] *= RAD2DEG;
+ }
+ td.tabscale *= RAD2DEG;
+ }
+ tab.n = td.nx;
+ tab.scale = td.tabscale;
+ snew(tab.data, tab.n*stride);
+ copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data);
+ done_tabledata(&td);
+
+ return tab;
+}
+
+t_forcetable *makeDispersionCorrectionTable(FILE *fp,
+ t_forcerec *fr, real rtab,
+ const char *tabfn)
+{
+ t_forcetable *dispersionCorrectionTable = NULL;
+
+ if (tabfn == NULL)
+ {
+ if (debug)
+ {
+ fprintf(debug, "No table file name passed, can not read table, can not do non-bonded interactions\n");
+ }
+ return dispersionCorrectionTable;
+ }
+
+ t_forcetable *fullTable = make_tables(fp, fr, tabfn, rtab, 0);
+ /* Copy the contents of the table to one that has just dispersion
+ * and repulsion, to improve cache performance. We want the table
+ * data to be aligned to 32-byte boundaries. The pointers could be
+ * freed but currently aren't. */
+ snew(dispersionCorrectionTable, 1);
+ dispersionCorrectionTable->interaction = GMX_TABLE_INTERACTION_VDWREP_VDWDISP;
+ dispersionCorrectionTable->format = fullTable->format;
+ dispersionCorrectionTable->r = fullTable->r;
+ dispersionCorrectionTable->n = fullTable->n;
+ dispersionCorrectionTable->scale = fullTable->scale;
+ dispersionCorrectionTable->formatsize = fullTable->formatsize;
+ dispersionCorrectionTable->ninteractions = 2;
+ dispersionCorrectionTable->stride = dispersionCorrectionTable->formatsize * dispersionCorrectionTable->ninteractions;
+ snew_aligned(dispersionCorrectionTable->data, dispersionCorrectionTable->stride*(dispersionCorrectionTable->n+1), 32);
+
+ for (int i = 0; i <= fullTable->n; i++)
+ {
+ for (int j = 0; j < 8; j++)
+ {
+ dispersionCorrectionTable->data[8*i+j] = fullTable->data[12*i+4+j];
+ }
+ }
+ sfree_aligned(fullTable->data);
+ sfree(fullTable);
+
+ return dispersionCorrectionTable;
+}
--- /dev/null
- bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle);
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2012,2014,2015,2016, 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.
+ */
+#ifndef GMX_TABLES_FORCETABLE_H
+#define GMX_TABLES_FORCETABLE_H
+
+#include <cstdio>
+
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/utility/real.h"
+
+#define GMX_MAKETABLES_FORCEUSER (1<<0)
+#define GMX_MAKETABLES_14ONLY (1<<1)
+
+/* Index in the tables that says which function to use */
+enum {
+ etiCOUL, etiLJ6, etiLJ12, etiNR
+};
+
+typedef double (*real_space_grid_contribution_computer)(double, double);
+/* Function pointer used to tell table_spline3_fill_ewald_lr whether it
+ * should calculate the grid contribution for electrostatics or LJ.
+ */
+
+void table_spline3_fill_ewald_lr(real *table_F,
+ real *table_V,
+ real *table_FDV0,
+ int ntab,
+ double dx,
+ real beta,
+ real_space_grid_contribution_computer v_lr);
+/* Fill tables of ntab points with spacing dr with the ewald long-range
+ * (mesh) force.
+ * There are three separate tables with format FDV0, F, and V.
+ * This function interpolates the Ewald mesh potential contribution
+ * with coefficient beta using a quadratic spline.
+ * The force can then be interpolated linearly.
+ */
+
+real ewald_spline3_table_scale(const interaction_const_t *ic);
+/* Return the scaling for the Ewald quadratic spline tables. */
+
+double v_q_ewald_lr(double beta, double r);
+/* Return the real space grid contribution for Ewald*/
+
+double v_lj_ewald_lr(double beta, double r);
+/* Return the real space grid contribution for LJ-Ewald*/
+
+t_forcetable *make_tables(FILE *fp,
+ const t_forcerec *fr,
+ const char *fn, real rtab, int flags);
+/* Return tables for inner loops. */
+
++bondedtable_t make_bonded_table(FILE *fplog, const char *fn, int angle);
+/* Return a table for bonded interactions,
+ * angle should be: bonds 0, angles 1, dihedrals 2
+ */
+
+/* Return a table for GB calculations */
+t_forcetable *make_gb_table(const t_forcerec *fr);
+
+/*! \brief Construct and return tabulated dispersion and repulsion interactions
+ *
+ * This table can be used to compute long-range dispersion corrections */
+t_forcetable *makeDispersionCorrectionTable(FILE *fp, t_forcerec *fr,
+ real rtab, const char *tabfn);
+
+#endif /* GMX_TABLES_FORCETABLE_H */
{
return str.find(substr) != std::string::npos;
}
+ //! \copydoc contains(const std::string &str, const char *substr)
+ static inline bool contains(const std::string &str, const std::string &substr)
+ {
+ return str.find(substr) != std::string::npos;
+ }
+/*!\brief Returns number of space-separated words in zero-terminated char ptr
+ *
+ * \param s Character pointer to zero-terminated, which will not be changed.
+ *
+ * \returns number of words in string.
+ *
+ * \note This routine is mainly meant to support legacy code in GROMACS. For
+ * new source you should try hard to use C++ string objects instead.
+ */
+std::size_t
+countWords(const char *s);
+
+/*!\brief Returns the number of space-separated words in a string object
+ *
+ * \param str Reference to string object, which will not be changed.
+ *
+ * \returns number of words in string.
+ */
+std::size_t
+countWords(const std::string &str);
+
++//! \copydoc endsWith(const std::string &str, const char *suffix)
++static inline bool endsWith(const std::string &str, const std::string &suffix)
++{
++ return endsWith(str, suffix.c_str());
++}
++
/*! \brief
* Removes a suffix from a string.
*
{ efXVG, "-dhdl", "dhdl", ffOPTWR },
{ efXVG, "-field", "field", ffOPTWR },
{ efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tabletf", "tabletf", ffOPTRD },
{ efXVG, "-tablep", "tablep", ffOPTRD },
- { efXVG, "-tableb", "table", ffOPTRD },
+ { efXVG, "-tableb", "table", ffOPTRDMULT },
{ efTRX, "-rerun", "rerun", ffOPTRD },
{ efXVG, "-tpi", "tpi", ffOPTWR },
{ efXVG, "-tpid", "tpidist", ffOPTWR },
fr = mk_forcerec();
fr->hwinfo = hwinfo;
fr->gpu_opt = &hw_opt->gpu_opt;
- init_forcerec(fplog, oenv, fr, fcd, inputrec, mtop, cr, box,
+ init_forcerec(fplog, fr, fcd, inputrec, mtop, cr, box,
opt2fn("-table", nfile, fnm),
- opt2fn("-tabletf", nfile, fnm),
opt2fn("-tablep", nfile, fnm),
- opt2fn("-tableb", nfile, fnm),
+ getFilenm("-tableb", nfile, fnm),
nbpu_opt,
FALSE,
pforce);
set(testname "MdrunTests")
set(exename "mdrun-test")
-gmx_build_unit_test(
- ${testname}
+gmx_add_gtest_executable(
${exename}
# files with code for tests
+ tabulated_bonded_interactions.cpp
+ energyreader.cpp
grompp.cpp
rerun.cpp
trajectory_writing.cpp
--- /dev/null
-#include "gromacs/utility/file.h"
+ /*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016, 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.
+ */
+
+ /*! \internal \file
+ * \brief
+ * Tests for tabulated bonded interactions
+ *
+ * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \ingroup module_mdrun_integration_tests
+ */
+ #include "gmxpre.h"
+
+ #include <gtest/gtest.h>
+
- File::writeFileFromString(runner_.topFileName_, formatString(g_butaneTopFileFormatString, interaction));
+ #include "gromacs/utility/stringutil.h"
++#include "gromacs/utility/textwriter.h"
+
+ #include "testutils/cmdlinetest.h"
+
+ #include "moduletest.h"
+
+ namespace gmx
+ {
+ namespace test
+ {
+
+ //! Format string for building a configurable .top file
+ const char *g_butaneTopFileFormatString = "\
+ [ defaults ]\n\
+ ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n\
+ 1 1 no 1.0 1.0\n\
+ \n\
+ [ atomtypes ]\n\
+ ;name mass charge ptype c6 c12\n\
+ CH2 14.02700 0.000 A 0.90975E-02 0.35333E-04\n\
+ CH3 15.03500 0.000 A 0.88765E-02 0.26150E-04\n\
+ \n\
+ [ moleculetype ]\n\
+ ; name nrexcl\n\
+ butane 3\n\
+ \n\
+ [ atoms ]\n\
+ ; nr type resnr residu atom cgnr\n\
+ 1 CH3 1 BUT C1 1\n\
+ 2 CH2 1 BUT C2 2\n\
+ 3 CH2 1 BUT C3 3\n\
+ 4 CH3 1 BUT C4 4\n\
+ \n\
+ %s\n\
+ \n\
+ [ system ]\n\
+ ; The name of the system to be simulated\n\
+ A single butane\n\
+ \n\
+ [ molecules ]\n\
+ ; Molname Number\n\
+ Butane 1\n\
+ ";
+
+ //! Test fixture for bonded interactions
+ class BondedInteractionsTest : public gmx::test::MdrunTestFixture
+ {
+ public:
+ //! Execute the trajectory writing test
+ void setupGrompp(const char *interaction)
+ {
+ runner_.topFileName_ = fileManager_.getTemporaryFilePath("butane1.top");
++ TextWriter::writeFileFromString(runner_.topFileName_, formatString(g_butaneTopFileFormatString, interaction));
+ runner_.groFileName_ = fileManager_.getInputFilePath("butane1.gro");
+ runner_.ndxFileName_ = fileManager_.getInputFilePath("butane1.ndx");
+ /* TODO Now that Verlet is the default, change the implementation
+ of useEmptyMdpFile() to do that. */
+ runner_.useStringAsMdpFile("");
+ }
+ //! Prepare an mdrun caller
+ CommandLine setupMdrun()
+ {
+ CommandLine rerunCaller;
+ rerunCaller.append("mdrun");
+ rerunCaller.addOption("-rerun", runner_.groFileName_);
+ return rerunCaller;
+ }
+ //! Check the output of mdrun
+ void checkMdrun()
+ {
+ // TODO verifying some energies and forces would be good,
+ // once other code in gerrit is reviewed
+ }
+ };
+
+ // This test ensures that a normal non-tabulated bond interaction works
+ TEST_F(BondedInteractionsTest, NormalBondWorks)
+ {
+ setupGrompp("[ bonds ]\n\
+ ; ai aj funct c0 c1\n\
+ 1 2 1 1.530000e-01 3.347000e+05");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ // This test ensures that a normal abulated bond interaction works
+ TEST_F(BondedInteractionsTest, TabulatedBondWorks)
+ {
+ setupGrompp("[ bonds ]\n\
+ ; ai aj funct n k\n\
+ 1 2 8 0 1000");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ std::string tableFileName = fileManager_.getInputFilePath("butane_b0.xvg");
+ rerunCaller.addOption("-tableb", tableFileName);
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ // This test ensures that a normal non-tabulated angle interaction works
+ TEST_F(BondedInteractionsTest, NormalAngleWorks)
+ {
+ setupGrompp("[ angles ]\n\
+ ; ai aj ak funct c0 c1\n\
+ 1 2 3 1 1.110000e+02 4.602000e+02");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ // This test ensures that a tabulated angle interaction works
+ TEST_F(BondedInteractionsTest, TabulatedAngleWorks)
+ {
+ setupGrompp("[ angles ]\n\
+ ; ai aj ak funct n k\n\
+ 1 2 3 8 0 1000");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ std::string tableFileName = fileManager_.getInputFilePath("butane_a0.xvg");
+ rerunCaller.addOption("-tableb", tableFileName);
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ // This test ensures that a normal non-tabulated dihedral interaction works
+ TEST_F(BondedInteractionsTest, NormalDihedralWorks)
+ {
+ setupGrompp("[ dihedrals ]\n \
+ ; ai aj ak al funct c0 c1 c2 c3 c4 c5\n\
+ 1 2 3 4 3 9.2789 12.156 -13.12 -3.0597 26.24 -31.495");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ // This test ensures that a tabulated dihedral interaction works
+ TEST_F(BondedInteractionsTest, TabulatedDihedralWorks)
+ {
+ setupGrompp("[ dihedrals ]\n\
+ ; ai aj ak al funct n k\n\
+ 1 2 3 4 8 0 1000");
+ EXPECT_EQ(0, runner_.callGrompp());
+
+ test::CommandLine rerunCaller = setupMdrun();
+ std::string tableFileName = fileManager_.getInputFilePath("butane_d0.xvg");
+ rerunCaller.addOption("-tableb", tableFileName);
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ checkMdrun();
+ }
+
+ } // namespace
+
+ } // namespace
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff