From: Mark Abraham Date: Tue, 7 Jun 2016 23:14:01 +0000 (+0200) Subject: Merge branch release-5-1 into release-2016 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=7139bb20270be53182fe4e6f34a4ca3c87159551;p=alexxy%2Fgromacs.git Merge branch release-5-1 into release-2016 Conflicts: CMakeLists.txt Adjacent unrelated changes, trivial resolution cmake/gmxDetectSimd.cmake Two unrelated fixes, managing the new GMX_STDLIB_LIBRARIES better, and working around the fact that we can't use try_run with noisy compilers. The new context invalidates the suggestion to use try_run once we require CMake 2.8.11. cmake/gmxSetBuildInformation.cmake As above. src/gromacs/commandline/filenm.h Adjacent changes caused by introducing Doxygen src/gromacs/fileio/filetypes.cpp As above src/gromacs/legacyheaders/force.h Change to function signature to pass const char *fn. src/gromacs/mdlib/forcerec.cpp Minor clashes from reorganized include files. TODO make_bonded_tables was refactored in release-5-1 to need t_filenm, which reintroduces a dependency on commandline/filenm.h. Decide what to do about this. src/gromacs/mdlib/mdatoms.cpp git was confused about the frozen-atom fix because of the indenting change, but the merge is straightforward in terms of code logic. src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda.cu New Doxygen in release-2016 clashed with introducing support for CUDA 6.0/6.1. src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh As above src/gromacs/tables/forcetable.cpp Change to function signature to pass const char *fn. src/programs/mdrun/tests/CMakeLists.txt Adjacent new files introduceed in both branches Change-Id: Iaaffacc186aa5ff67c83522d2c07b05afeec75b2 --- 7139bb20270be53182fe4e6f34a4ca3c87159551 diff --cc CMakeLists.txt index b05d6867f2,4fee4802b6..849383540c --- a/CMakeLists.txt +++ b/CMakeLists.txt @@@ -564,11 -489,12 +565,15 @@@ if (TMPI_ATOMICS_DISABLED add_definitions(-DTMPI_ATOMICS_DISABLED) endif() +# Note this relies on zlib detection having already run +include(gmxManageTNG) + - # now that we have detected the dependencies, do the second GPU configure pass - gmx_gpu_setup() + 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) diff --cc cmake/gmxDetectSimd.cmake index 4abf65e290,78225473f9..28fbdc21f6 --- a/cmake/gmxDetectSimd.cmake +++ b/cmake/gmxDetectSimd.cmake @@@ -68,75 -72,55 +72,57 @@@ function(gmx_suggest_simd _suggested_si 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) diff --cc cmake/gmxSetBuildInformation.cmake index bf6551f9e8,c53f369d5b..e3d677ac57 --- a/cmake/gmxSetBuildInformation.cmake +++ b/cmake/gmxSetBuildInformation.cmake @@@ -81,89 -81,108 +81,108 @@@ macro(gmx_set_build_information 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) diff --cc src/gromacs/commandline/filenm.cpp index b7f7d8b5a0,0000000000..f8b5e55282 mode 100644,000000..100644 --- a/src/gromacs/commandline/filenm.cpp +++ b/src/gromacs/commandline/filenm.cpp @@@ -1,291 -1,0 +1,306 @@@ +/* + * 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, by the GROMACS development team, led by ++ * 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 +#include + +#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; + } +} diff --cc src/gromacs/commandline/filenm.h index 19c8f744de,0000000000..1ca8470663 mode 100644,000000..100644 --- a/src/gromacs/commandline/filenm.h +++ b/src/gromacs/commandline/filenm.h @@@ -1,179 -1,0 +1,185 @@@ +/* + * 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, by the GROMACS development team, led by ++ * 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 diff --cc src/gromacs/fileio/filetypes.cpp index 9886c35b56,0000000000..2002498f34 mode 100644,000000..100644 --- a/src/gromacs/fileio/filetypes.cpp +++ b/src/gromacs/fileio/filetypes.cpp @@@ -1,323 -1,0 +1,323 @@@ +/* + * 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, by the GROMACS development team, led by ++ * 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 + +#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; +} diff --cc src/gromacs/mdlib/forcerec.cpp index 6ed1608b20,d5710ab84b..703a1cd471 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@@ -44,23 -42,32 +44,24 @@@ #include #include +#include + #include ++#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" @@@ -85,13 -80,9 +86,14 @@@ #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" @@@ -2265,8 -2358,9 +2321,8 @@@ void init_forcerec(FILE *f 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) @@@ -3049,20 -3213,26 +3105,26 @@@ 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 { diff --cc src/gromacs/mdlib/forcerec.h index 9ff5b250c1,0000000000..274d909bf8 mode 100644,000000..100644 --- a/src/gromacs/mdlib/forcerec.h +++ b/src/gromacs/mdlib/forcerec.h @@@ -1,137 -1,0 +1,138 @@@ +/* + * 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, by the GROMACS development team, led by ++ * 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 - * \param[in] tabbfn Table potential file for bonded interactions ++ * \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 char *tabbfn, ++ 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 diff --cc src/gromacs/mdlib/mdatoms.cpp index ffa5825bec,0000000000..63ff5e8b49 mode 100644,000000..100644 --- a/src/gromacs/mdlib/mdatoms.cpp +++ b/src/gromacs/mdlib/mdatoms.cpp @@@ -1,435 -1,0 +1,438 @@@ +/* + * 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) 2012,2013,2014,2015, by the GROMACS development team, led by ++ * 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 + +#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. - * Note that constraints can still move a partially frozen particle. ++ /* 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; +} diff --cc src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda.cu index 3946419231,0429445dd0..058680e11d --- a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda.cu +++ b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda.cu @@@ -71,11 -78,15 +71,12 @@@ texture 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; diff --cc src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh index c63eda7f43,91ab0c6967..de26e947b7 --- a/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh +++ b/src/gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel.cuh @@@ -109,26 -109,22 +109,28 @@@ * 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) diff --cc src/gromacs/tables/forcetable.cpp index fa6a198749,0000000000..2b52c178a3 mode 100644,000000..100644 --- a/src/gromacs/tables/forcetable.cpp +++ b/src/gromacs/tables/forcetable.cpp @@@ -1,1611 -1,0 +1,1614 @@@ +/* + * 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 + +#include + +#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(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; - ssd += fabs(2*(f - numf)/(f + numf)); ++ if (f + numf != 0) ++ { ++ ssd += fabs(2*(f - numf)/(f + numf)); ++ } + ns++; + } + } + if (ns > 0) + { + ssd /= ns; - 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)); ++ 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 rrc, 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(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(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, char *fn, int angle) ++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; +} diff --cc src/gromacs/tables/forcetable.h index 7e68fdaf2c,0000000000..9236b34405 mode 100644,000000..100644 --- a/src/gromacs/tables/forcetable.h +++ b/src/gromacs/tables/forcetable.h @@@ -1,101 -1,0 +1,101 @@@ +/* + * 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 + +#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, char *fn, int angle); ++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 */ diff --cc src/gromacs/utility/stringutil.h index 864f58f8b9,b2b790e502..f54aedaa9e --- a/src/gromacs/utility/stringutil.h +++ b/src/gromacs/utility/stringutil.h @@@ -114,28 -115,12 +114,39 @@@ static inline bool contains(const std:: { 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. * diff --cc src/programs/mdrun/mdrun.cpp index 9458dfbf91,fcdc36952a..e56c4a1453 --- a/src/programs/mdrun/mdrun.cpp +++ b/src/programs/mdrun/mdrun.cpp @@@ -238,8 -238,9 +239,8 @@@ int gmx_mdrun(int argc, char *argv[] { 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 }, diff --cc src/programs/mdrun/runner.cpp index 3294776ff1,7489b07dcb..969a3d9513 --- a/src/programs/mdrun/runner.cpp +++ b/src/programs/mdrun/runner.cpp @@@ -1150,10 -1097,11 +1150,10 @@@ int mdrunner(gmx_hw_opt_t *hw_opt 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); diff --cc src/programs/mdrun/tests/CMakeLists.txt index c84a5d074d,22dc2532fc..cfbdb09568 --- a/src/programs/mdrun/tests/CMakeLists.txt +++ b/src/programs/mdrun/tests/CMakeLists.txt @@@ -35,10 -35,11 +35,11 @@@ 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 diff --cc src/programs/mdrun/tests/tabulated_bonded_interactions.cpp index 0000000000,9a7f331507..4e18b667f8 mode 000000,100644..100644 --- a/src/programs/mdrun/tests/tabulated_bonded_interactions.cpp +++ b/src/programs/mdrun/tests/tabulated_bonded_interactions.cpp @@@ -1,0 -1,209 +1,209 @@@ + /* + * 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 + * \ingroup module_mdrun_integration_tests + */ + #include "gmxpre.h" + + #include + -#include "gromacs/utility/file.h" + #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"); - File::writeFileFromString(runner_.topFileName_, formatString(g_butaneTopFileFormatString, interaction)); ++ 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