# - Check the username performing the build, as well as date and time
#
-# gmx_detect_simd(GMX_SUGGESTED_SIMD)
+# gmx_detect_simd(_suggested_simd_)
#
-# Try to detect CPU information and suggest SIMD instruction set
+# Try to detect CPU features and suggest a SIMD instruction set
# that fits the current CPU. This should work on all architectures
# where we are not cross-compiling; depending on the architecture the
# detection will either use special assembly instructions (like cpuid),
# preprocessor defines, or probing /proc/cpuinfo on Linux.
#
-# This assumes gmx_detect_target_architecture() has already been run,
-# so that things like GMX_TARGET_X86 are already available.
-# (otherwise we cannot use inline ASM on x86).
-#
-# Sets ${GMX_SUGGESTED_SIMD} in the parent scope if
-# GMX_SIMD is not set (e.g. by the user, or a previous run
-# of CMake).
+# Sets ${suggested_simd} in the parent scope if GMX_SIMD is not set
+# (e.g. by the user, or a previous run of CMake).
#
# we rely on inline asm support for GNU!
include(gmxTestInlineASM)
+# Ensure things like GMX_TARGET_X86 are available
+include(gmxDetectTargetArchitecture)
+gmx_detect_target_architecture()
+include(gmxDetectCpu)
function(gmx_suggest_simd _suggested_simd)
- if(${_suggested_simd})
- # There's already been a suggestion made, which can't change
- return()
- endif()
-
- # for x86 we need inline asm to use cpuid
- gmx_test_inline_asm_gcc_x86(GMX_X86_GCC_INLINE_ASM)
-
- if(GMX_X86_GCC_INLINE_ASM)
- set(GCC_INLINE_ASM_DEFINE "-DGMX_X86_GCC_INLINE_ASM=1")
- else()
- set(GCC_INLINE_ASM_DEFINE "-DGMX_X86_GCC_INLINE_ASM=0")
+ if (NOT SUGGEST_SIMD_QUIETLY)
+ message(STATUS "Detecting best SIMD instructions for this CPU")
endif()
- message(STATUS "Detecting best SIMD instructions for this CPU")
-
- # Get CPU SIMD properties information
- if(GMX_TARGET_X86)
- set(GMX_TARGET_X86_VALUE 1)
- else()
- set(GMX_TARGET_X86_VALUE 0)
- endif()
- set(_compile_definitions "${GCC_INLINE_ASM_DEFINE} -I${CMAKE_SOURCE_DIR}/src -DGMX_CPUINFO_STANDALONE ${GMX_STDLIB_CXX_FLAGS} -DGMX_TARGET_X86=${GMX_TARGET_X86_VALUE}")
-
# 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)
- # TODO Extract this try_compile to a helper function, because
- # it duplicates code in gmxSetBuildInformation.cmake
- set(GMX_DETECTSIMD_BINARY "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/GmxDetectSimd${CMAKE_EXECUTABLE_SUFFIX}")
- set(LINK_LIBRARIES "${GMX_STDLIB_LIBRARIES}")
- try_compile(GMX_DETECTSIMD_COMPILED
- "${CMAKE_CURRENT_BINARY_DIR}"
- "${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)
+ # Detect CPU features and place the string in CPU_DETECTION_FEATURES
+ # Note that we are NOT limited to x86.
+ gmx_run_cpu_detection(features)
+
+ if (DEFINED CPU_DETECTION_FEATURES)
+ # Make a concrete suggestion of SIMD level if a feature flag
+ # matches. Make sure that the match strings below work even if
+ # the feature is first or last.
+ set(CPU_DETECTION_FEATURES " ${CPU_DETECTION_FEATURES} ")
- if(GMX_DETECTSIMD_COMPILED)
- if(NOT DEFINED GMX_DETECTSIMD_RUN)
- execute_process(COMMAND ${GMX_DETECTSIMD_BINARY} "-features"
- RESULT_VARIABLE GMX_DETECTSIMD_RUN
- OUTPUT_VARIABLE OUTPUT_TMP
- ERROR_QUIET)
- set(GMX_DETECTSIMD_RUN "${GMX_DETECTSIMD_RUN}" CACHE INTERNAL "Result of running cpuinfo code to detect SIMD support")
- if(GMX_DETECTSIMD_RUN EQUAL 0)
- # Make a concrete suggestion of SIMD level
- 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")
- else()
- # Intel
- set(OUTPUT_SIMD "AVX_256")
- 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 " AND NOT GMX_DOUBLE)
- set(OUTPUT_SIMD "ARM_NEON")
- endif()
- endif()
- message(STATUS "Detected best SIMD instructions for this CPU - ${OUTPUT_SIMD}")
+ if(GMX_TARGET_X86)
+ if(CPU_DETECTION_FEATURES MATCHES " avx512er ")
+ set(OUTPUT_SIMD "AVX_512_KNL")
+ elseif(CPU_DETECTION_FEATURES MATCHES " avx512f ")
+ set(OUTPUT_SIMD "AVX_512")
+ elseif(CPU_DETECTION_FEATURES MATCHES " avx2 ")
+ if(CPU_DETECTION_FEATURES MATCHES " amd ")
+ set(OUTPUT_SIMD "AVX2_128")
+ else()
+ set(OUTPUT_SIMD "AVX2_256")
+ endif()
+ elseif(CPU_DETECTION_FEATURES MATCHES " avx ")
+ if(CPU_DETECTION_FEATURES MATCHES " fma4 ")
+ # AMD that works better with avx-128-fma
+ set(OUTPUT_SIMD "AVX_128_FMA")
else()
- message(WARNING "Cannot run cpuinfo code, which means no SIMD suggestion can be made.")
- message(STATUS "Run output: ${OUTPUT_TMP}")
+ # Intel
+ set(OUTPUT_SIMD "AVX_256")
endif()
+ elseif(CPU_DETECTION_FEATURES MATCHES " sse4.1 ")
+ set(OUTPUT_SIMD "SSE4.1")
+ elseif(CPU_DETECTION_FEATURES MATCHES " sse2 ")
+ set(OUTPUT_SIMD "SSE2")
endif()
else()
- message(WARNING "Cannot compile cpuinfo code, which means no SIMD instructions.")
- message(STATUS "Compile output: ${GMX_DETECTSIMD_COMPILED_OUTPUT}")
+ if(CPU_DETECTION_FEATURES MATCHES " vsx ")
+ set(OUTPUT_SIMD "IBM_VSX")
+ elseif(CPU_DETECTION_FEATURES MATCHES " vmx ")
+ set(OUTPUT_SIMD "IBM_VMX")
+ elseif(CPU_DETECTION_FEATURES MATCHES " qpx ")
+ set(OUTPUT_SIMD "IBM_QPX")
+ elseif(CPU_DETECTION_FEATURES MATCHES " neon_asimd ")
+ set(OUTPUT_SIMD "ARM_NEON_ASIMD")
- elseif(CPU_DETECTION_FEATURES MATCHES " neon ")
++ elseif(CPU_DETECTION_FEATURES MATCHES " neon " AND NOT GMX_DOUBLE)
+ set(OUTPUT_SIMD "ARM_NEON")
+ endif()
+ endif()
+ if (NOT SUGGEST_SIMD_QUIETLY)
+ message(STATUS "Detected best SIMD instructions for this CPU - ${OUTPUT_SIMD}")
endif()
else()
- message(WARNING "Cannot detect SIMD architecture for this cross-compile; you should check it manually.")
+ if (NOT SUGGEST_SIMD_QUIETLY)
+ message(STATUS "Detection for best SIMD instructions failed, using SIMD - ${OUTPUT_SIMD}")
+ endif()
endif()
- set(${_suggested_simd} "${OUTPUT_SIMD}" CACHE INTERNAL "Suggested SIMD")
+ set(${_suggested_simd} "${OUTPUT_SIMD}" PARENT_SCOPE)
+ set(SUGGEST_SIMD_QUIETLY TRUE CACHE INTERNAL "Be quiet during future construction of SIMD suggestions")
endfunction()
function(gmx_detect_simd _suggested_simd)
- if(NOT DEFINED GMX_SIMD)
+ if(GMX_SIMD STREQUAL "AUTO")
if(GMX_TARGET_BGQ)
+ # BG/Q requires cross-compilation, so needs this
+ # logic. While the qpx feature flag in cpuinfo works, it
+ # can't be returned by cpuinfo running on the build host.
set(${_suggested_simd} "IBM_QPX")
elseif(GMX_TARGET_FUJITSU_SPARC64)
# HPC-ACE is always present. In the future we
# pain as much as possible:
# - use the CUDA_HOST_COMPILER if defined by the user, otherwise
# - auto-detect compatible nvcc host compiler and set nvcc -ccbin (if not MPI wrapper)
-# - set icc compatibility mode to gcc 4.6
+# - set icc compatibility mode to gcc 4.8.1
# - (advanced) variables set:
# * CUDA_HOST_COMPILER - the host compiler for nvcc (only with cmake <2.8.10)
# * CUDA_HOST_COMPILER_OPTIONS - the full host-compiler related option list passed to nvcc
unset(CUDA_HOST_COMPILER_AUTOSET CACHE)
endif()
-# Set the host compiler for nvcc if this is not set by CMake (v<=2.8.9)
-#
-# Note that even though nvcc compiles host code as C++, we use the
-# CMAKE_C_COMPILER as host compiler. We do this because CUDA versions
-# preceding 5.0 only recognize icc, but not icpc. However, both gcc and icc
-# (i.e. all supported compilers) happily compile C++ code.
-#
-# Also note that with MSVC nvcc sets the -compiler-bindir option behind the
-# scenes; to avoid conflicts we don't set -ccbin automatically.
-#
-# TODO: remove this when CMAke >=v2.8.10 is required.
-if (NOT DEFINED CUDA_HOST_COMPILER AND NOT MSVC)
- set(CUDA_HOST_COMPILER "${CMAKE_C_COMPILER}")
- set(CUDA_HOST_COMPILER_AUTOSET TRUE CACHE INTERNAL
- "True if CUDA_HOST_COMPILER is automatically set")
-endif()
-
# glibc 2.23 changed string.h in a way that breaks CUDA compilation in
# many projects, but which has a trivial workaround. It would be nicer
# to compile with nvcc and see that the workaround is necessary and
# set up host compiler and its options
if(CUDA_HOST_COMPILER_CHANGED)
- # FindCUDA in CMake 2.8.10 sets the host compiler internally
- if (CMAKE_VERSION VERSION_LESS "2.8.10")
- set(CUDA_HOST_COMPILER ${CUDA_HOST_COMPILER}
- CACHE PATH "Host compiler for nvcc")
- endif()
-
- # On *nix force icc in gcc 4.6 compatibility mode. This is needed
- # as even with icc used as host compiler, when icc's gcc compatibility
- # mode is higher than the max gcc version officially supported by CUDA,
- # nvcc will freak out.
set(CUDA_HOST_COMPILER_OPTIONS "")
- if (UNIX AND
- ((CMAKE_C_COMPILER_ID MATCHES "Intel" AND
- (CUDA_HOST_COMPILER_AUTOSET OR CMAKE_C_COMPILER STREQUAL CUDA_HOST_COMPILER)) OR
- (CMAKE_CXX_COMPILER_ID MATCHES "Intel" AND CMAKE_CXX_COMPILER STREQUAL CUDA_HOST_COMPILER))
- )
- message(STATUS "Setting Intel Compiler compatibity mode to gcc 4.6 for nvcc host compilation")
- list(APPEND CUDA_HOST_COMPILER_OPTIONS "-Xcompiler;-gcc-version=460")
- endif()
if(APPLE AND CMAKE_C_COMPILER_ID MATCHES "GNU")
# Some versions of gcc-4.8 and gcc-4.9 produce errors (in particular on OS X)
# => compile sm_20, sm_30, sm_35, sm_37, sm_50, & sm_52 SASS, and compute_52 PTX
# - with CUDA >=8.0 CC 6.0-6.2 is supported (but we know nothing about CC 6.2, so we won't generate code or it)
# => compile sm_20, sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_61 SASS, and compute_60 and compute_61 PTX
- #
+ # - with CUDA >=9.0 CC 7.0 is supported and CC 2.0 is no longer supported
+ # => compile sm_30, sm_35, sm_37, sm_50, sm_52, sm_60, sm_61, sm_70 SASS, and compute_70 PTX
#
# Note that CUDA 6.5.19 second patch release supports cc 5.2 too, but
# CUDA_VERSION does not contain patch version and having PTX 5.0 JIT-ed is
# equally fast as compiling with sm_5.2 anyway.
# First add flags that trigger SASS (binary) code generation for physical arch
- list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_20,code=sm_20")
+ if(CUDA_VERSION VERSION_LESS "9.00") # < 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_20,code=sm_20")
+ endif()
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_30,code=sm_30")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_35,code=sm_35")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_60,code=sm_60")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_61,code=sm_61")
endif()
+ if(NOT CUDA_VERSION VERSION_LESS "9.0") # >= 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_70,code=sm_70")
+ endif()
# Next add flags that trigger PTX code generation for the newest supported virtual arch
# that's useful to JIT to future architectures
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_50,code=compute_50")
elseif(CUDA_VERSION VERSION_LESS "8.0")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_52,code=compute_52")
- else() # version >= 8.0
+ elseif(CUDA_VERSION VERSION_LESS "9.0")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_60,code=compute_60")
list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_61,code=compute_61")
+ else() # version >= 9.0
+ list (APPEND GMX_CUDA_NVCC_GENCODE_FLAGS "-gencode;arch=compute_70,code=compute_70")
endif()
endif()
-gmx_dependent_cache_variable(GMX_CUDA_TARGET_SM "List of CUDA GPU architecture codes to compile for (without the sm_ prefix)" STRING "" GMX_CUDA_TARGET_SM)
-gmx_dependent_cache_variable(GMX_CUDA_TARGET_COMPUTE "List of CUDA virtual architecture codes to compile for (without the compute_ prefix)" STRING "" GMX_CUDA_TARGET_COMPUTE)
+if (GMX_CUDA_TARGET_SM)
+ set_property(CACHE GMX_CUDA_TARGET_SM PROPERTY HELPSTRING "List of CUDA GPU architecture codes to compile for (without the sm_ prefix)")
+ set_property(CACHE GMX_CUDA_TARGET_SM PROPERTY TYPE STRING)
+endif()
+if (GMX_CUDA_TARGET_COMPUTE)
+ set_property(CACHE GMX_CUDA_TARGET_COMPUTE PROPERTY HELPSTRING "List of CUDA virtual architecture codes to compile for (without the compute_ prefix)")
+ set_property(CACHE GMX_CUDA_TARGET_COMPUTE PROPERTY TYPE STRING)
+endif()
# assemble the CUDA flags
list(APPEND GMX_CUDA_NVCC_FLAGS "${GMX_CUDA_NVCC_GENCODE_FLAGS}")
endif()
# assemble the CUDA host compiler flags
-# with CMake <2.8.10 the host compiler needs to be set on the nvcc command line
-if (CMAKE_VERSION VERSION_LESS "2.8.10")
- list(APPEND GMX_CUDA_NVCC_FLAGS "-ccbin=${CUDA_HOST_COMPILER}")
-endif()
list(APPEND GMX_CUDA_NVCC_FLAGS "${CUDA_HOST_COMPILER_OPTIONS}")
# The flags are set as local variables which shadow the cache variables. The cache variables
# (can be set by the user) are appended. This is done in a macro to set the flags when all
# host compiler flags are already set.
macro(GMX_SET_CUDA_NVCC_FLAGS)
- if(CUDA_PROPAGATE_HOST_FLAGS)
- set(CUDA_PROPAGATE_HOST_FLAGS OFF)
-
- # When CUDA 6.5 is required we should use C++11 also for CUDA and also propagate
- # the C++11 flag to CUDA. Then we can use the solution implemented in FindCUDA
- # (starting with 3.3 - can be backported). For now we need to remove the C++11
- # flag which means we need to manually propagate all other flags.
- string(REGEX REPLACE "[-]+std=c\\+\\+0x" "" _CMAKE_CXX_FLAGS_SANITIZED "${CMAKE_CXX_FLAGS}")
-
- # The IBM xlc compiler chokes if we use both altivec and Cuda. Solve
- # this by not propagating the flag in this case.
- if(CMAKE_CXX_COMPILER_ID MATCHES "XL")
- string(REGEX REPLACE "-qaltivec" "" _CMAKE_CXX_FLAGS_SANITIZED "${_CMAKE_CXX_FLAGS_SANITIZED}")
- endif()
-
- # CUDA versions prior to 7.5 come with a header (math_functions.h) which uses the _MSC_VER macro
- # unconditionally, so we strip -Wundef from the propagatest flags for earlier CUDA versions.
- if (CUDA_VERSION VERSION_LESS "7.5")
- string(REGEX REPLACE "-Wundef" "" _CMAKE_CXX_FLAGS_SANITIZED "${_CMAKE_CXX_FLAGS_SANITIZED}")
- endif()
-
- string(REPLACE " " "," _flags "${_CMAKE_CXX_FLAGS_SANITIZED}")
- set(CUDA_NVCC_FLAGS "${GMX_CUDA_NVCC_FLAGS};${CUDA_NVCC_FLAGS};-Xcompiler;${_flags}")
-
- # Create list of all possible configurations. For multi-configuration this is CMAKE_CONFIGURATION_TYPES
- # and for single configuration CMAKE_BUILD_TYPE. Not sure why to add the default ones, but FindCUDA
- # claims one should.
- set(CUDA_configuration_types ${CMAKE_CONFIGURATION_TYPES} ${CMAKE_BUILD_TYPE} Debug MinSizeRel Release RelWithDebInfo)
- list(REMOVE_DUPLICATES CUDA_configuration_types)
-
- foreach(_config ${CUDA_configuration_types})
- string(TOUPPER ${_config} _config_upper)
- string(REPLACE " " "," _flags "${CMAKE_CXX_FLAGS_${_config_upper}}")
- set(CUDA_NVCC_FLAGS_${_config_upper} "${CUDA_NVCC_FLAGS_${_config_upper}};-Xcompiler;${_flags}")
- endforeach()
- else()
- set(CUDA_NVCC_FLAGS "${GMX_CUDA_NVCC_FLAGS};${CUDA_NVCC_FLAGS}")
- endif()
+ set(CUDA_NVCC_FLAGS "${GMX_CUDA_NVCC_FLAGS};${CUDA_NVCC_FLAGS}")
endmacro()
# 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 avx test source, used if the AVX flags are set below
-include(gmxTestAVXMaskload)
+include(gmxDetectCpu)
include(gmxFindFlagsForSource)
# Macro that manages setting the respective C and C++ toolchain
# VSX uses the same function API as Altivec/VMX, so make sure we tune for the current CPU and not VMX.
# By putting these flags here rather than in the general compiler flags file we can safely assume
# that we are at least on Power7 since that is when VSX appeared.
- if(BUILD_CPU_BRAND MATCHES "POWER7")
+ gmx_run_cpu_detection(brand)
+ if(CPU_DETECTION_BRAND MATCHES "POWER7")
gmx_test_cflag(GNU_C_VSX_POWER7 "-mcpu=power7 -mtune=power7" ${TOOLCHAIN_C_FLAGS_VARIABLE})
gmx_test_cflag(GNU_CXX_VSX_POWER7 "-mcpu=power7 -mtune=power7" ${TOOLCHAIN_CXX_FLAGS_VARIABLE})
else()
#
# Section to set (and test) compiler flags for SIMD.
#
-# The flags will be set based on the GMX_SIMD choice provided by the user.
-# Automatic detection of the architecture on the build host is done prior to
-# calling this macro.
+# If the user chose the (default) automatic behaviour, then detection
+# is run to suggest a SIMD choice suitable for the build
+# host. Otherwise, the users's choice is always honoured. The compiler
+# flags will be set based on that choice.
#
-if(GMX_SIMD STREQUAL "NONE")
+set(GMX_SIMD_ACTIVE ${GMX_SIMD})
+if(GMX_SIMD STREQUAL "AUTO")
+ include(gmxDetectSimd)
+ gmx_detect_simd(GMX_SUGGESTED_SIMD)
+ set(GMX_SIMD_ACTIVE ${GMX_SUGGESTED_SIMD})
+endif()
+
+if(GMX_SIMD_ACTIVE STREQUAL "NONE")
# nothing to do configuration-wise
set(SIMD_STATUS_MESSAGE "SIMD instructions disabled")
-elseif(GMX_SIMD STREQUAL "SSE2")
+elseif(GMX_SIMD_ACTIVE STREQUAL "SSE2")
gmx_find_flags(
"#include<xmmintrin.h>
int main(){__m128 x=_mm_set1_ps(0.5);x=_mm_rsqrt_ps(x);return _mm_movemask_ps(x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-msse2" "/arch:SSE2" "-hgnu")
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("SSE2" "disable SIMD support (slow)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling SSE2 SIMD instructions")
-elseif(GMX_SIMD STREQUAL "SSE4.1")
+elseif(GMX_SIMD_ACTIVE STREQUAL "SSE4.1")
# Note: MSVC enables SSE4.1 with the SSE2 flag, so we include that in testing.
gmx_find_flags(
set(GMX_SIMD_X86_SSE4_1 1)
set(SIMD_STATUS_MESSAGE "Enabling SSE4.1 SIMD instructions")
-elseif(GMX_SIMD STREQUAL "AVX_128_FMA")
+elseif(GMX_SIMD_ACTIVE STREQUAL "AVX_128_FMA")
prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
${INCLUDE_INTRIN_H}
int main(){__m128 x=_mm_set1_ps(0.5);x=_mm_macc_ps(x,x,x);return _mm_movemask_ps(x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-mfma4" "-hgnu")
# We only need to check the last (FMA) test; that will always fail if the generic AVX test failed
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("128-bit AVX with FMA support" "choose SSE4.1 SIMD (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling 128-bit AVX SIMD GROMACS SIMD (with fused-multiply add)")
- gmx_test_avx_gcc_maskload_bug(GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG "${SIMD_C_FLAGS}")
-
-elseif(GMX_SIMD STREQUAL "AVX_256")
+elseif(GMX_SIMD_ACTIVE STREQUAL "AVX_256")
prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
"#include<immintrin.h>
int main(){__m256 x=_mm256_set1_ps(0.5);x=_mm256_add_ps(x,x);return _mm256_movemask_ps(x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-mavx" "/arch:AVX" "-hgnu")
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("AVX" "choose SSE4.1 SIMD (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling 256-bit AVX SIMD instructions")
- gmx_test_avx_gcc_maskload_bug(GMX_SIMD_X86_AVX_GCC_MASKLOAD_BUG "${SIMD_C_FLAGS}")
-
-elseif(GMX_SIMD STREQUAL "AVX2_256")
+elseif(GMX_SIMD_ACTIVE MATCHES "AVX2_")
prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
"#include<immintrin.h>
int main(){__m256i x=_mm256_set1_epi32(5);x=_mm256_add_epi32(x,x);return _mm256_movemask_epi8(x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-march=core-avx2" "-mavx2" "/arch:AVX" "-hgnu") # no AVX2-specific flag for MSVC yet
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("AVX2" "choose AVX SIMD (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
- set(SIMD_STATUS_MESSAGE "Enabling 256-bit AVX2 SIMD instructions")
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
- # No need to test for Maskload bug - it was fixed before gcc added AVX2 support
+ if(GMX_SIMD_ACTIVE STREQUAL "AVX2_128")
+ set(SIMD_STATUS_MESSAGE "Enabling 128-bit AVX2 SIMD instructions")
+ else()
+ set(SIMD_STATUS_MESSAGE "Enabling 256-bit AVX2 SIMD instructions")
+ endif()
-elseif(GMX_SIMD STREQUAL "MIC")
+elseif(GMX_SIMD_ACTIVE STREQUAL "MIC")
# No flags needed. Not testing.
set(GMX_SIMD_X86_MIC 1)
set(SIMD_STATUS_MESSAGE "Enabling MIC (Xeon Phi) SIMD instructions")
-elseif(GMX_SIMD STREQUAL "AVX_512")
+elseif(GMX_SIMD_ACTIVE STREQUAL "AVX_512")
prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
"#include<immintrin.h>
int main(){__m512 y,x=_mm512_set1_ps(0.5);y=_mm512_fmadd_ps(x,x,x);return (int)_mm512_cmp_ps_mask(x,y,_CMP_LT_OS);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-xCORE-AVX512" "-mavx512f -mfma" "-mavx512f" "/arch:AVX" "-hgnu") # no AVX_512F flags known for MSVC yet
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("AVX 512F" "choose a lower level of SIMD (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling 512-bit AVX-512 SIMD instructions")
-elseif(GMX_SIMD STREQUAL "AVX_512_KNL")
+elseif(GMX_SIMD_ACTIVE STREQUAL "AVX_512_KNL")
prepare_x86_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
"#include<immintrin.h>
int main(){__m512 y,x=_mm512_set1_ps(0.5);y=_mm512_rsqrt28_ps(x);return (int)_mm512_cmp_ps_mask(x,y,_CMP_LT_OS);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-xMIC-AVX512" "-mavx512er -mfma" "-mavx512er" "/arch:AVX" "-hgnu") # no AVX_512ER flags known for MSVC yet
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("AVX 512ER" "choose a lower level of SIMD (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_X86_${GMX_SIMD} 1)
+ set(GMX_SIMD_X86_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling 512-bit AVX-512-KNL SIMD instructions")
-elseif(GMX_SIMD STREQUAL "ARM_NEON")
+elseif(GMX_SIMD_ACTIVE STREQUAL "ARM_NEON")
+ if (GMX_DOUBLE)
+ message(FATAL_ERROR "ARM_NEON SIMD support is not available for a double precision build because the architecture lacks double-precision support")
+ endif()
+
gmx_find_flags(
"#include<arm_neon.h>
int main(){float32x4_t x=vdupq_n_f32(0.5);x=vmlaq_f32(x,x,x);return vgetq_lane_f32(x,0)>0;}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-mfpu=neon-vfpv4" "-mfpu=neon" "")
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("ARM NEON" "disable SIMD support (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling 32-bit ARM NEON SIMD instructions")
-elseif(GMX_SIMD STREQUAL "ARM_NEON_ASIMD")
+elseif(GMX_SIMD_ACTIVE STREQUAL "ARM_NEON_ASIMD")
gmx_find_flags(
"#include<arm_neon.h>
int main(){float64x2_t x=vdupq_n_f64(0.5);x=vfmaq_f64(x,x,x);x=vrndnq_f64(x);return vgetq_lane_f64(x,0)>0;}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"")
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("ARM (AArch64) NEON Advanced SIMD" "particularly gcc version 4.9 or later, or disable SIMD support (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling ARM (AArch64) NEON Advanced SIMD instructions")
-elseif(GMX_SIMD STREQUAL "IBM_QPX")
+elseif(GMX_SIMD_ACTIVE STREQUAL "IBM_QPX")
try_compile(TEST_QPX ${CMAKE_BINARY_DIR}
"${CMAKE_SOURCE_DIR}/cmake/TestQPX.c")
if (TEST_QPX)
message(WARNING "IBM QPX SIMD instructions selected. This will work, but SIMD kernels are only available for the Verlet cut-off scheme. The plain C kernels that are used for the group cut-off scheme kernels will be slow, so please consider using the Verlet cut-off scheme.")
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling IBM QPX SIMD instructions")
else()
- gmx_give_fatal_error_when_simd_support_not_found("IBM QPX" "or 'cmake .. -DCMAKE_TOOLCHAIN_FILE=Platform/BlueGeneQ-static-XL-CXX' to set up the tool chain" "${SUGGEST_BINUTILS_UPDATE}")
+ gmx_give_fatal_error_when_simd_support_not_found("IBM QPX" "or 'cmake .. -DCMAKE_TOOLCHAIN_FILE=Platform/BlueGeneQ-static-bgclang-CXX' to set up the tool chain" "${SUGGEST_BINUTILS_UPDATE}")
endif()
-elseif(GMX_SIMD STREQUAL "IBM_VMX")
+elseif(GMX_SIMD_ACTIVE STREQUAL "IBM_VMX")
gmx_find_flags(
"#include<altivec.h>
int main(){vector float x,y=vec_ctf(vec_splat_s32(1),0);x=vec_madd(y,y,y);return vec_all_ge(y,x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-maltivec -mabi=altivec" "-qarch=auto -qaltivec")
- if(NOT SIMD_${GMX_SIMD}_C_FLAGS OR NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS OR NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("IBM VMX" "disable SIMD support (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling IBM VMX SIMD instructions")
-elseif(GMX_SIMD STREQUAL "IBM_VSX")
+elseif(GMX_SIMD_ACTIVE STREQUAL "IBM_VSX")
prepare_power_vsx_toolchain(TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS)
"#include<altivec.h>
int main(){vector double x,y=vec_splats(1.0);x=vec_madd(y,y,y);return vec_all_ge(y,x);}"
TOOLCHAIN_C_FLAGS TOOLCHAIN_CXX_FLAGS
- SIMD_${GMX_SIMD}_C_FLAGS SIMD_${GMX_SIMD}_CXX_FLAGS
+ SIMD_${GMX_SIMD_ACTIVE}_C_FLAGS SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS
"-mvsx" "-maltivec -mabi=altivec" "-qarch=auto -qaltivec")
# Usually we check also for the C compiler here, but a C compiler
# at least version 3.7 cannot pass this check with the C compiler
# in the latest xlc 13.1.5, but the C++ compiler has different
# behaviour and is OK. See Redmine #2102.
- if(NOT SIMD_${GMX_SIMD}_CXX_FLAGS)
+ if(NOT SIMD_${GMX_SIMD_ACTIVE}_CXX_FLAGS)
gmx_give_fatal_error_when_simd_support_not_found("IBM VSX" "disable SIMD support (slower)" "${SUGGEST_BINUTILS_UPDATE}")
endif()
set(SIMD_C_FLAGS "${TOOLCHAIN_C_FLAGS}")
set(SIMD_CXX_FLAGS "${TOOLCHAIN_CXX_FLAGS}")
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling IBM VSX SIMD instructions")
-elseif(GMX_SIMD STREQUAL "SPARC64_HPC_ACE")
+elseif(GMX_SIMD_ACTIVE STREQUAL "SPARC64_HPC_ACE")
# Note that GMX_RELAXED_DOUBLE_PRECISION is enabled by default in the top-level CMakeLists.txt
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling Sparc64 HPC-ACE SIMD instructions")
-elseif(GMX_SIMD STREQUAL "REFERENCE")
+elseif(GMX_SIMD_ACTIVE STREQUAL "REFERENCE")
# NB: This file handles settings for the SIMD module, so in the interest
# of proper modularization, please do NOT put any verlet kernel settings in this file.
add_definitions(-DGMX_SIMD_REF_DOUBLE_WIDTH=${GMX_SIMD_REF_DOUBLE_WIDTH})
endif()
- set(GMX_SIMD_${GMX_SIMD} 1)
+ set(GMX_SIMD_${GMX_SIMD_ACTIVE} 1)
set(SIMD_STATUS_MESSAGE "Enabling reference (emulated) SIMD instructions.")
else()
- gmx_invalid_option_value(GMX_SIMD)
+ gmx_invalid_option_value(GMX_SIMD_ACTIVE)
endif()
-gmx_check_if_changed(SIMD_CHANGED GMX_SIMD)
+gmx_check_if_changed(SIMD_CHANGED GMX_SIMD_ACTIVE)
if (SIMD_CHANGED AND DEFINED SIMD_STATUS_MESSAGE)
message(STATUS "${SIMD_STATUS_MESSAGE}")
endif()
if(NOT DEFINED GMX_SIMD_CALLING_CONVENTION)
if(GMX_TARGET_BGQ)
set(CALLCONV_LIST " ")
- elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND GMX_SIMD STREQUAL "REFERENCE")
+ elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND GMX_SIMD_ACTIVE STREQUAL "REFERENCE")
set(CALLCONV_LIST __regcall " ")
elseif(CMAKE_CXX_COMPILER_ID MATCHES "XL")
set(CALLCONV_LIST " ")
endforeach()
endif()
+if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+ # GCC bug 49001, 54412 on Windows (just warn, since it might be fixed in later versions)
+ if((CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.9.0" OR CMAKE_SIZEOF_VOID_P EQUAL 8)
+ AND (WIN32 OR CYGWIN)
+ AND (GMX_SIMD_ACTIVE MATCHES "AVX") AND NOT (GMX_SIMD_ACTIVE STREQUAL "AVX_128_FMA"))
+ message(WARNING "GCC on Windows (GCC older than 4.9 in 32-bit mode, or any version in 64-bit mode) with 256-bit AVX will probably crash. You might want to choose a different GMX_SIMD or a different compiler.")
+ endif()
+endif()
+
endmacro()
user-guide/mdrun-features.rst
user-guide/mdrun-performance.rst
user-guide/mdp-options.rst
+ user-guide/run-time-errors.rst
user-guide/file-formats.rst
user-guide/cmdline.rst
user-guide/environment-variables.rst
+ user-guide/terminology.rst
user-guide/plotje.gif
user-guide/xvgr.gif
conf.py
EXTRA_VARS
SPHINX_EXTENSION_PATH RELENG_PATH
EXPECTED_DOXYGEN_VERSION
- GMX_CMAKE_MINIMUM_REQUIRED_VERSION REQUIRED_CUDA_VERSION
+ CMAKE_MINIMUM_REQUIRED_VERSION REQUIRED_CUDA_VERSION
REQUIRED_OPENCL_MIN_VERSION
REQUIRED_CUDA_COMPUTE_CAPABILITY REGRESSIONTEST_VERSION
SOURCE_MD5SUM REGRESSIONTEST_MD5SUM_STRING
comprehensive comparisons have been carried out, and some caution must
be taking in interpreting the results.
+When using temperature and/or pressure coupling the total energy is
+no longer conserved. Instead there is a \normindex{conserved energy quantity}
+the formula of which will depend on the combination or temperature and
+pressure coupling algorithm used. For all coupling algorithms, except
+for Andersen temperature coupling and Parrinello-Rahman pressure coupling
+combined with shear stress, the conserved energy quantity is computed
+and stored in the energy and log file. Note that this quantity will not
+be conserved when external forces are applied to the system, such as
+pulling on group with a changing distance or an electric field.
+Furthermore, how well the energy is conserved depends on the accuracy
+of all algorithms involved in the simulation. Usually the algorithms that
+cause most drift are constraints and the pair-list buffer, depending
+on the parameters used.
+
\subsubsection{Berendsen temperature coupling\pawsindexquiet{Berendsen}{temperature coupling}\index{weak coupling}}
The Berendsen algorithm mimics weak coupling with first-order
kinetics to an external heat bath with given temperature $T_0$.
numbers which may crash the simulation. In normal use,
$\lambda$ will always be much closer to 1.0.
+The thermostat modifies the kinetic energy at each scaling step by:
+\beq
+\Delta E_k = (\lambda - 1)^2 E_k
+\eeq
+The sum of these changes over the run needs to subtracted from the total energy
+to obtain the conserved energy quantity.
+
\subsubsection{Velocity-rescaling temperature coupling\pawsindexquiet{velocity-rescaling}{temperature coupling}}
The velocity-rescaling thermostat~\cite{Bussi2007a} is essentially a Berendsen
thermostat (see above) with an additional stochastic term that ensures
This thermostat produces a correct canonical ensemble and still has
the advantage of the Berendsen thermostat: first order decay of
temperature deviations and no oscillations.
-When an $NVT$ ensemble is used, the conserved energy quantity
-is written to the energy and log file.
\subsubsection{\normindex{Andersen thermostat}}
One simple way to maintain a thermostatted ensemble is to take an
\end{array}\right).
\eeq
The velocities are neither scaled nor rotated.
+Since the equations of motion are modified by pressure coupling, the conserved
+energy quantity also needs to be modified. For first order pressure coupling,
+the work the barostat applies to the system every step needs to
+be subtracted from the total energy to obtain the conserved energy quantity:
+\beq
+- \sum_{i,j} (\mu_{ij} -\delta_{ij}) P_{ij} V =
+\sum_{i,j} 2(\mu_{ij} -\delta_{ij}) \Xi_{ij}
+\eeq
+where $\delta_{ij}$ is the Kronecker delta and ${\bf \Xi}$ is the virial.
+Note that the factor 2 originates from the factor $\frac{1}{2}$
+in the virial definition (\eqnref{Xi}).
+
In {\gromacs}, the Berendsen scaling can also be done isotropically,
which means that instead of $\ve{P}$ a diagonal matrix with elements of size
for the Nos{\'e}-Hoover coupling. In most cases you would combine the
Parrinello-Rahman barostat with the Nos{\'e}-Hoover
thermostat, but to keep it simple we only show the Parrinello-Rahman
-modification here:
-
+modification here. The modified Hamiltonian, which will be conserved, is:
+\beq
+E_\mathrm{pot} + E_\mathrm{kin} + \sum_i P_{ii} V +
+\sum_{i,j} \frac{1}{2} W_{ij} \left( \frac{\de b_{ij}}{\de t} \right)^2
+\eeq
+The equations of motion for the atoms, obtained from the Hamiltonian are:
\bea \frac {\de^2\ve{r}_i}{\de t^2} & = & \frac{\ve{F}_i}{m_i} -
\ve{M} \frac{\de \ve{r}_i}{\de t} , \\ \ve{M} & = & \ve{b}^{-1} \left[
\ve{b} \frac{\de \ve{b}'}{\de t} + \frac{\de \ve{b}}{\de t} \ve{b}'
the run input file, one should use the minimized conformation from
the full precision trajectory file, as the structure file is not
accurate enough.
- {\tt \normindex{g_nmeig}} does the diagonalization and
+ {\tt \normindex{gmx nmeig}} does the diagonalization and
the sorting of the normal modes according to their frequencies.
- Both {\tt mdrun} and {\tt g_nmeig} should be run in double precision.
- The normal modes can be analyzed with the program {\tt g_anaeig}.
+ Both {\tt mdrun} and {\tt gmx nmeig} should be run in double precision.
+ The normal modes can be analyzed with the program {\tt gmx anaeig}.
Ensembles of structures at any temperature and for any subset of
- normal modes can be generated with {\tt \normindex{g_nmens}}.
+ normal modes can be generated with {\tt \normindex{gmx nmens}}.
An overview of normal-mode analysis and the related principal component
analysis (see \secref{covanal}) can be found in~\cite{Hayward95b}.
% } % Brace matches ifthenelse test for gmxlite
{\gromacs} now also supports the use of Bennett's Acceptance Ratio~\cite{Bennett1976}
for calculating values of $\Delta$G for transformations from state A to state B using
- the program {\tt \normindex{g_bar}}. The same data can also be used to calculate free
+ the program {\tt \normindex{gmx bar}}. The same data can also be used to calculate free
energies using MBAR~\cite{Shirts2008}, though the analysis currently requires external tools from
the external {\tt pymbar} package, at https://SimTK.org/home/pymbar.
The procedure for essential dynamics sampling or flooding is as follows.
First, the eigenvectors and eigenvalues need to be determined
- using covariance analysis ({\tt g_covar})
- or normal-mode analysis ({\tt g_nmeig}).
+ using covariance analysis ({\tt gmx covar})
+ or normal-mode analysis ({\tt gmx nmeig}).
Then, this information is fed into {\tt make_edi},
which has many options for selecting vectors and setting parameters,
see {\tt gmx make_edi -h}.
%
% This file is part of the GROMACS molecular simulation package.
%
- % Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ % Copyright (c) 2013,2014,2015,2016,2017, 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.
local rotation axis of the slab and the Gaussian-weighted positions.
+\section{\normindex{Electric fields}}
+A pulsed and oscillating electric field can be applied according to:
+\begin{equation}
+E(t) = E_0 \exp\left[-\frac{(t-t_0)^2}{2\sigma^2}\right]\cos\left[\omega (t-t_0)\right]
+\label{eq_efield}
+\end{equation}
+where $E_0$ is the field strength, the angular frequency \mbox{$\omega = 2\pi c/\lambda$}, $t_0$ is
+the time at of the peak in the field strength and $\sigma$ is the with
+of the pulse. Special cases occur when $\sigma$ = 0 (non-pulsed field)
+and for $\omega$ is 0 (static field).
+
+This simulated \normindex{laser}-pulse was applied to
+simulations of melting ice~\cite{Caleman2008a}. A pulsed electric field may
+look ike Fig.~\ref{fig:field}. In the supporting
+information of that paper the impact of an applied electric field on a
+system under periodic boundary conditions is analyzed. It is described
+that the effective electric field under PBC is larger than the applied
+field, by a factor depending on the size of the box and the dielectric
+properties of molecules in the box. For a system with static dielectric
+properties this factor can be corrected for. But for a system where
+the dielectric varies over time, for example a membrane protein with
+a pore that opens and closes during the simulatippn, this way of applying
+an electric field is not useful. In such cases one can use the computational
+electrophysiology protocol described in the next section (\secref{compel}).
+\begin{figure}[ht]
+\centerline{\includegraphics[width=8cm]{plots/field}}
+\caption {A simulated laser pulse in GROMACS.}
+\label{fig:field}
+\end{figure}
+
+Electric fields are applied when the following options are specified
+in the {\tt grompp.mdp} file. You specify, in order, $E_0$, $\omega$,
+$t_0$ and $\sigma$:
+\begin{verbatim}
+ElectricField-x = 0.04 0 0 0
+\end{verbatim}
+yields a static field with $E_0$ = 0.04 V/nm in the X-direction. In contrast,
+\begin{verbatim}
+ElectricField-x = 2.0 150 5 0
+\end{verbatim}
+yields an oscillating electric field with $E_0$ = 2 V/nm, $\omega$ = 150/ps and
+$t_0$ = 5 ps. Finally
+\begin{verbatim}
+ElectricField-x = 2.0 150 5 1
+\end{verbatim}
+yields an pulsed-oscillating electric field with $E_0$ = 2 V/nm, $\omega$ = 150/ps and
+$t_0$ = 5 ps and $\sigma$ = 1 ps. Read more in ref.~\cite{Caleman2008a}.
+Note that the input file format is changed from the undocumented older
+version. A figure like Fig.~\ref{fig:field} may be produced by passing
+the {\tt -field} option to {\tt gmx mdrun}.
+
+
\section{\normindex{Computational Electrophysiology}}
\label{sec:compel}
\left( \int_{t_0}^{{t_0}+t} P_{xz}(t') \mbox{d} t' \right)^2
\right\rangle_{t_0}
\eeq
- This can be done with {\tt g_energy}.
+ This can be done with {\tt gmx energy}.
This method converges very slowly~\cite{Hess2002a}, and as such
a nanosecond simulation might not
be long enough for an accurate determination of the viscosity.
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/ewald/pme.h"
#include "gromacs/fft/calcgrid.h"
-#include "gromacs/gmxlib/md_logging.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
+#include "gromacs/mdlib/nbnxn_pairlist.h"
#include "gromacs/mdlib/sim_util.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/smalloc.h"
#include "pme-internal.h"
/*! \brief Parameters and settings for one PP-PME setup */
struct pme_setup_t {
real rcut_coulomb; /**< Coulomb cut-off */
- real rlist; /**< pair-list cut-off */
+ real rlistOuter; /**< cut-off for the outer pair-list */
+ real rlistInner; /**< cut-off for the inner pair-list */
real spacing; /**< (largest) PME grid spacing */
ivec grid; /**< the PME grid dimensions */
real grid_efficiency; /**< ineffiency factor for non-uniform grids <= 1 */
real cut_spacing; /**< the minimum cutoff / PME grid spacing ratio */
real rcut_vdw; /**< Vdw cutoff (does not change) */
real rcut_coulomb_start; /**< Initial electrostatics cutoff */
- real rbuf_coulomb; /**< the pairlist buffer size */
- real rbuf_vdw; /**< the pairlist buffer size */
+ real rbufOuter_coulomb; /**< the outer pairlist buffer size */
+ real rbufOuter_vdw; /**< the outer pairlist buffer size */
+ real rbufInner_coulomb; /**< the inner pairlist buffer size */
+ real rbufInner_vdw; /**< the inner pairlist buffer size */
matrix box_start; /**< the initial simulation box */
int n; /**< the count of setup as well as the allocation size */
pme_setup_t *setup; /**< the PME+cutoff setups */
* read bActive anywhere */
bool pme_loadbal_is_active(const pme_load_balancing_t *pme_lb)
{
- return pme_lb != NULL && pme_lb->bActive;
+ return pme_lb != nullptr && pme_lb->bActive;
}
void pme_loadbal_init(pme_load_balancing_t **pme_lb_p,
t_commrec *cr,
- FILE *fp_log,
+ const gmx::MDLogger &mdlog,
const t_inputrec *ir,
matrix box,
const interaction_const_t *ic,
- struct gmx_pme_t *pmedata,
+ const NbnxnListParameters *listParams,
+ gmx_pme_t *pmedata,
gmx_bool bUseGPU,
gmx_bool *bPrinting)
{
+ GMX_RELEASE_ASSERT(ir->cutoff_scheme != ecutsGROUP, "PME tuning is not supported with cutoff-scheme=group (because it contains bugs)");
+
pme_load_balancing_t *pme_lb;
real spm, sp;
int d;
snew(pme_lb, 1);
- pme_lb->bSepPMERanks = !(cr->duty & DUTY_PME);
+ pme_lb->bSepPMERanks = !(cr->duty & DUTY_PME);
/* Initially we turn on balancing directly on based on PP/PME imbalance */
- pme_lb->bTriggerOnDLB = FALSE;
+ pme_lb->bTriggerOnDLB = FALSE;
/* Any number of stages >= 2 is supported */
- pme_lb->nstage = 2;
+ pme_lb->nstage = 2;
- pme_lb->cutoff_scheme = ir->cutoff_scheme;
+ pme_lb->cutoff_scheme = ir->cutoff_scheme;
- pme_lb->rbuf_coulomb = ic->rlist - ic->rcoulomb;
- pme_lb->rbuf_vdw = ic->rlist - ic->rvdw;
+ pme_lb->rbufOuter_coulomb = listParams->rlistOuter - ic->rcoulomb;
+ pme_lb->rbufOuter_vdw = listParams->rlistOuter - ic->rvdw;
+ pme_lb->rbufInner_coulomb = listParams->rlistInner - ic->rcoulomb;
+ pme_lb->rbufInner_vdw = listParams->rlistInner - ic->rvdw;
copy_mat(box, pme_lb->box_start);
if (ir->ePBC == epbcXY && ir->nwall == 2)
pme_lb->cur = 0;
pme_lb->setup[0].rcut_coulomb = ic->rcoulomb;
- pme_lb->setup[0].rlist = ic->rlist;
+ pme_lb->setup[0].rlistOuter = listParams->rlistOuter;
+ pme_lb->setup[0].rlistInner = listParams->rlistInner;
pme_lb->setup[0].grid[XX] = ir->nkx;
pme_lb->setup[0].grid[YY] = ir->nky;
pme_lb->setup[0].grid[ZZ] = ir->nkz;
if (!wallcycle_have_counter())
{
- md_print_warn(cr, fp_log, "NOTE: Cycle counters unsupported or not enabled in kernel. Cannot use PME-PP balancing.\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("NOTE: Cycle counters unsupported or not enabled in kernel. Cannot use PME-PP balancing.");
}
/* Tune with GPUs and/or separate PME ranks.
dd_dlb_lock(cr->dd);
if (dd_dlb_is_locked(cr->dd))
{
- md_print_warn(cr, fp_log, "NOTE: DLB will not turn on during the first phase of PME tuning\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("NOTE: DLB will not turn on during the first phase of PME tuning");
}
}
real fac, sp;
real tmpr_coulomb, tmpr_vdw;
int d;
- gmx_bool grid_ok;
+ bool grid_ok;
/* Try to add a new setup with next larger cut-off to the list */
pme_lb->n++;
srenew(pme_lb->setup, pme_lb->n);
set = &pme_lb->setup[pme_lb->n-1];
- set->pmedata = NULL;
+ set->pmedata = nullptr;
get_pme_nnodes(dd, &npmeranks_x, &npmeranks_y);
fac *= 1.01;
clear_ivec(set->grid);
- sp = calc_grid(NULL, pme_lb->box_start,
- fac*pme_lb->setup[pme_lb->cur].spacing,
- &set->grid[XX],
- &set->grid[YY],
- &set->grid[ZZ]);
+ sp = calcFftGrid(nullptr, pme_lb->box_start,
+ fac*pme_lb->setup[pme_lb->cur].spacing,
+ minimalPmeGridSize(pme_order),
+ &set->grid[XX],
+ &set->grid[YY],
+ &set->grid[ZZ]);
/* As here we can't easily check if one of the PME ranks
* uses threading, we do a conservative grid check.
* This means we can't use pme_order or less grid lines
* per PME rank along x, which is not a strong restriction.
*/
- gmx_pme_check_restrictions(pme_order,
- set->grid[XX], set->grid[YY], set->grid[ZZ],
- npmeranks_x, npmeranks_y,
- TRUE,
- FALSE,
- &grid_ok);
+ grid_ok = gmx_pme_check_restrictions(pme_order,
+ set->grid[XX], set->grid[YY], set->grid[ZZ],
+ npmeranks_x,
+ true,
+ false);
}
while (sp <= 1.001*pme_lb->setup[pme_lb->cur].spacing || !grid_ok);
if (pme_lb->cutoff_scheme == ecutsVERLET)
{
/* Never decrease the Coulomb and VdW list buffers */
- set->rlist = std::max(set->rcut_coulomb + pme_lb->rbuf_coulomb,
- pme_lb->rcut_vdw + pme_lb->rbuf_vdw);
+ set->rlistOuter = std::max(set->rcut_coulomb + pme_lb->rbufOuter_coulomb,
+ pme_lb->rcut_vdw + pme_lb->rbufOuter_vdw);
+ set->rlistInner = std::max(set->rcut_coulomb + pme_lb->rbufInner_coulomb,
+ pme_lb->rcut_vdw + pme_lb->rbufInner_vdw);
}
else
{
- tmpr_coulomb = set->rcut_coulomb + pme_lb->rbuf_coulomb;
- tmpr_vdw = pme_lb->rcut_vdw + pme_lb->rbuf_vdw;
+ /* TODO Remove these lines and pme_lb->cutoff_scheme */
- set->rlist = std::min(tmpr_coulomb, tmpr_vdw);
+ tmpr_coulomb = set->rcut_coulomb + pme_lb->rbufOuter_coulomb;
+ tmpr_vdw = pme_lb->rcut_vdw + pme_lb->rbufOuter_vdw;
+ /* Two (known) bugs with cutoff-scheme=group here:
+ * - This modification of rlist results in incorrect DD comunication.
+ * - We should set fr->bTwinRange = (fr->rlistlong > fr->rlist).
+ */
+ set->rlistOuter = std::min(tmpr_coulomb, tmpr_vdw);
+ set->rlistInner = set->rlistOuter;
}
- set->spacing = sp;
+ set->spacing = sp;
/* The grid efficiency is the size wrt a grid with uniform x/y/z spacing */
set->grid_efficiency = 1;
for (d = 0; d < DIM; d++)
pre,
desc, set->grid[XX], set->grid[YY], set->grid[ZZ], set->rcut_coulomb,
buft);
- if (fp_err != NULL)
+ if (fp_err != nullptr)
{
fprintf(fp_err, "\r%s\n", buf);
fflush(fp_err);
}
- if (fp_log != NULL)
+ if (fp_log != nullptr)
{
fprintf(fp_log, "%s\n", buf);
}
gmx_step_str(step, sbuf),
pmelblim_str[pme_lb->elimited],
pme_lb->setup[pme_loadbal_end(pme_lb)-1].rcut_coulomb);
- if (fp_err != NULL)
+ if (fp_err != nullptr)
{
fprintf(fp_err, "\r%s\n", buf);
fflush(fp_err);
}
- if (fp_log != NULL)
+ if (fp_log != nullptr)
{
fprintf(fp_log, "%s\n", buf);
}
t_commrec *cr,
FILE *fp_err,
FILE *fp_log,
+ const gmx::MDLogger &mdlog,
const t_inputrec *ir,
t_state *state,
double cycles,
* better overal performance can be obtained with a slightly
* shorter cut-off and better DD load balancing.
*/
- set_dd_dlb_max_cutoff(cr, pme_lb->setup[pme_lb->fastest].rlist);
+ set_dd_dlb_max_cutoff(cr, pme_lb->setup[pme_lb->fastest].rlistOuter);
}
}
cycles_fast = pme_lb->setup[pme_lb->fastest].cycles;
if (OK && ir->ePBC != epbcNONE)
{
- OK = (gmx::square(pme_lb->setup[pme_lb->cur+1].rlist)
+ OK = (gmx::square(pme_lb->setup[pme_lb->cur+1].rlistOuter)
<= max_cutoff2(ir->ePBC, state->box));
if (!OK)
{
if (DOMAINDECOMP(cr))
{
OK = change_dd_cutoff(cr, state, ir,
- pme_lb->setup[pme_lb->cur].rlist);
+ pme_lb->setup[pme_lb->cur].rlistOuter);
if (!OK)
{
/* Failed: do not use this setup */
if (DOMAINDECOMP(cr) && pme_lb->stage > 0)
{
- OK = change_dd_cutoff(cr, state, ir, pme_lb->setup[pme_lb->cur].rlist);
+ OK = change_dd_cutoff(cr, state, ir, pme_lb->setup[pme_lb->cur].rlistOuter);
if (!OK)
{
/* For some reason the chosen cut-off is incompatible with DD.
/* This should not happen, as we set limits on the DLB bounds.
* But we implement a complete failsafe solution anyhow.
*/
- md_print_warn(cr, fp_log, "The fastest PP/PME load balancing setting (cutoff %.3f nm) is no longer available due to DD DLB or box size limitations\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+ "The fastest PP/PME load balancing setting (cutoff %.3f nm) is no longer available due to DD DLB or box size limitations", pme_lb->fastest);
pme_lb->fastest = pme_lb->lower_limit;
pme_lb->start = pme_lb->lower_limit;
}
set = &pme_lb->setup[pme_lb->cur];
- ic->rcoulomb = set->rcut_coulomb;
- ic->rlist = set->rlist;
- ic->ewaldcoeff_q = set->ewaldcoeff_q;
+ NbnxnListParameters *listParams = nbv->listParams.get();
+
+ ic->rcoulomb = set->rcut_coulomb;
+ listParams->rlistOuter = set->rlistOuter;
+ listParams->rlistInner = set->rlistInner;
+ ic->ewaldcoeff_q = set->ewaldcoeff_q;
/* TODO: centralize the code that sets the potentials shifts */
if (ic->coulomb_modifier == eintmodPOTSHIFT)
{
- ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+ GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
+ ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb;
}
if (EVDW_PME(ic->vdwtype))
{
}
/* We always re-initialize the tables whether they are used or not */
- init_interaction_const_tables(NULL, ic, rtab);
+ init_interaction_const_tables(nullptr, ic, rtab);
- nbnxn_gpu_pme_loadbal_update_param(nbv, ic);
+ nbnxn_gpu_pme_loadbal_update_param(nbv, ic, listParams);
/* With tMPI + GPUs some ranks may be sharing GPU(s) and therefore
* also sharing texture references. To keep the code simple, we don't
if (!pme_lb->bSepPMERanks)
{
- if (pme_lb->setup[pme_lb->cur].pmedata == NULL)
+ if (pme_lb->setup[pme_lb->cur].pmedata == nullptr)
{
/* Generate a new PME data structure,
* copying part of the old pointers.
*/
gmx_pme_reinit(&set->pmedata,
cr, pme_lb->setup[0].pmedata, ir,
- set->grid);
+ set->grid, set->ewaldcoeff_q, set->ewaldcoeff_lj);
}
*pmedata = set->pmedata;
}
if (debug)
{
- print_grid(NULL, debug, "", "switched to", set, -1);
+ print_grid(nullptr, debug, "", "switched to", set, -1);
}
if (pme_lb->stage == pme_lb->nstage)
t_commrec *cr,
FILE *fp_err,
FILE *fp_log,
+ const gmx::MDLogger &mdlog,
const t_inputrec *ir,
t_forcerec *fr,
t_state *state,
{
/* Unlock the DLB=auto, DLB is allowed to activate */
dd_dlb_unlock(cr->dd);
- md_print_warn(cr, fp_log, "NOTE: DLB can now turn on, when beneficial\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("NOTE: DLB can now turn on, when beneficial");
/* We don't deactivate the tuning yet, since we will balance again
* after DLB gets turned on, if it does within PMETune_period.
* This also ensures that we won't disable the currently
* optimal setting during a second round of PME balancing.
*/
- set_dd_dlb_max_cutoff(cr, fr->ic->rlist);
+ set_dd_dlb_max_cutoff(cr, fr->nbv->listParams->rlistOuter);
}
}
* but the first data collected is skipped anyhow.
*/
pme_load_balance(pme_lb, cr,
- fp_err, fp_log,
+ fp_err, fp_log, mdlog,
ir, state, pme_lb->cycles_c - cycles_prev,
fr->ic, fr->nbv, &fr->pmedata,
step);
/* Update constants in forcerec/inputrec to keep them in sync with fr->ic */
fr->ewaldcoeff_q = fr->ic->ewaldcoeff_q;
fr->ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
- fr->rlist = fr->ic->rlist;
+ fr->rlist = fr->nbv->listParams->rlistOuter;
fr->rcoulomb = fr->ic->rcoulomb;
fr->rvdw = fr->ic->rvdw;
if (ir->eDispCorr != edispcNO)
{
- calc_enervirdiff(NULL, ir->eDispCorr, fr);
+ calc_enervirdiff(nullptr, ir->eDispCorr, fr);
}
}
{
/* Make sure DLB is allowed when we deactivate PME tuning */
dd_dlb_unlock(cr->dd);
- md_print_warn(cr, fp_log, "NOTE: DLB can now turn on, when beneficial\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("NOTE: DLB can now turn on, when beneficial");
}
*bPrinting = pme_lb->bBalance;
fprintf(fplog,
" %-7s %6.3f nm %6.3f nm %3d %3d %3d %5.3f nm %5.3f nm\n",
name,
- setup->rcut_coulomb, setup->rlist,
+ setup->rcut_coulomb, setup->rlistInner,
setup->grid[XX], setup->grid[YY], setup->grid[ZZ],
setup->spacing, 1/setup->ewaldcoeff_q);
}
/*! \brief Print all load-balancing settings */
static void print_pme_loadbal_settings(pme_load_balancing_t *pme_lb,
- t_commrec *cr,
FILE *fplog,
+ const gmx::MDLogger &mdlog,
gmx_bool bNonBondedOnGPU)
{
double pp_ratio, grid_ratio;
real pp_ratio_temporary;
- pp_ratio_temporary = pme_lb->setup[pme_lb->cur].rlist / pme_lb->setup[0].rlist;
+ pp_ratio_temporary = pme_lb->setup[pme_lb->cur].rlistInner / pme_lb->setup[0].rlistInner;
pp_ratio = gmx::power3(pp_ratio_temporary);
grid_ratio = pme_grid_points(&pme_lb->setup[pme_lb->cur])/
(double)pme_grid_points(&pme_lb->setup[0]);
if (pp_ratio > 1.5 && !bNonBondedOnGPU)
{
- md_print_warn(cr, fplog,
- "NOTE: PME load balancing increased the non-bonded workload by more than 50%%.\n"
- " For better performance, use (more) PME ranks (mdrun -npme),\n"
- " or if you are beyond the scaling limit, use fewer total ranks (or nodes).\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText(
+ "NOTE: PME load balancing increased the non-bonded workload by more than 50%.\n"
+ " For better performance, use (more) PME ranks (mdrun -npme),\n"
+ " or if you are beyond the scaling limit, use fewer total ranks (or nodes).");
}
else
{
}
void pme_loadbal_done(pme_load_balancing_t *pme_lb,
- t_commrec *cr,
FILE *fplog,
+ const gmx::MDLogger &mdlog,
gmx_bool bNonBondedOnGPU)
{
- if (fplog != NULL && (pme_lb->cur > 0 || pme_lb->elimited != epmelblimNO))
+ if (fplog != nullptr && (pme_lb->cur > 0 || pme_lb->elimited != epmelblimNO))
{
- print_pme_loadbal_settings(pme_lb, cr, fplog, bNonBondedOnGPU);
+ print_pme_loadbal_settings(pme_lb, fplog, mdlog, bNonBondedOnGPU);
}
/* TODO: Here we should free all pointers in pme_lb,
curr->delta_t = dt;
curr->beginfit = (1 - 2*GMX_REAL_EPS)*beginfit;
curr->endfit = (1 + 2*GMX_REAL_EPS)*endfit;
- curr->x0 = NULL;
- curr->n_offs = NULL;
+ curr->x0 = nullptr;
+ curr->n_offs = nullptr;
curr->nframes = 0;
curr->nlast = 0;
curr->dim_factor = dim_factor;
}
for (i = 0; (i < nrgrp); i++)
{
- curr->ndata[i] = NULL;
- curr->data[i] = NULL;
+ curr->ndata[i] = nullptr;
+ curr->data[i] = nullptr;
if (bTen)
{
- curr->datam[i] = NULL;
+ curr->datam[i] = nullptr;
}
}
- curr->time = NULL;
- curr->lsq = NULL;
+ curr->time = nullptr;
+ curr->lsq = nullptr;
curr->nmol = nmol;
if (curr->nmol > 0)
{
if (DD)
{
fprintf(out, "# MSD gathered over %g %s with %d restarts\n",
- msdtime, output_env_get_time_unit(oenv), curr->nrestart);
+ msdtime, output_env_get_time_unit(oenv).c_str(), curr->nrestart);
fprintf(out, "# Diffusion constants fitted from time %g to %g %s\n",
- beginfit, endfit, output_env_get_time_unit(oenv));
+ beginfit, endfit, output_env_get_time_unit(oenv).c_str());
for (i = 0; i < curr->ngrp; i++)
{
fprintf(out, "# D[%10s] = %.4f (+/- %.4f) (1e-5 cm^2/s)\n",
gmx_stats_t lsq1;
int i, j;
real a, b, D, Dav, D2av, VarD, sqrtD, sqrtD_max, scale;
- t_pdbinfo *pdbinfo = NULL;
- const int *mol2a = NULL;
+ t_pdbinfo *pdbinfo = nullptr;
+ const int *mol2a = nullptr;
out = xvgropen(fn, "Diffusion Coefficients / Molecule", "Molecule", "D", oenv);
gmx_stats_add_point(lsq1, xx, yy, dx, dy);
}
}
- gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, NULL, NULL, NULL, NULL);
+ gmx_stats_get_ab(lsq1, elsqWEIGHT_NONE, &a, &b, nullptr, nullptr, nullptr, nullptr);
gmx_stats_free(lsq1);
D = a*FACTOR/curr->dim_factor;
if (D < 0)
{
scale *= 10;
}
- GMX_RELEASE_ASSERT(pdbinfo != NULL, "Internal error - pdbinfo not set for PDB input");
+ GMX_RELEASE_ASSERT(pdbinfo != nullptr, "Internal error - pdbinfo not set for PDB input");
for (i = 0; i < top->atoms.nr; i++)
{
pdbinfo[i].bfac *= scale;
}
- write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, NULL, ePBC, box);
+ write_sto_conf(fn_pdb, "molecular MSD", &top->atoms, x, nullptr, ePBC, box);
}
}
#define prev (1-cur)
matrix box;
gmx_bool bFirst;
- gmx_rmpbc_t gpbc = NULL;
+ gmx_rmpbc_t gpbc = nullptr;
natoms = read_first_x(oenv, &status, fn, &curr->t0, &(x[cur]), box);
#ifdef DEBUG
fprintf(stderr, "Read %d atoms for first frame\n", natoms);
#endif
- if ((gnx_com != NULL) && natoms < top->atoms.nr)
+ if ((gnx_com != nullptr) && natoms < top->atoms.nr)
{
fprintf(stderr, "WARNING: The trajectory only contains part of the system (%d of %d atoms) and therefore the COM motion of only this part of the system will be removed\n", natoms, top->atoms.nr);
}
t = curr->t0;
if (x_pdb)
{
- *x_pdb = NULL;
+ *x_pdb = nullptr;
}
if (bMol)
t_pdb > t - 0.5*(t - t_prev) &&
t_pdb < t + 0.5*(t - t_prev))))
{
- if (*x_pdb == NULL)
+ if (*x_pdb == nullptr)
{
snew(*x_pdb, natoms);
}
{
for (i = 0; (i < curr->ngrp); i++)
{
- curr->ndata[i] = NULL;
- curr->data[i] = NULL;
+ curr->ndata[i] = nullptr;
+ curr->data[i] = nullptr;
if (bTen)
{
- curr->datam[i] = NULL;
+ curr->datam[i] = nullptr;
}
}
- curr->time = NULL;
+ curr->time = nullptr;
}
maxframes += 10;
for (i = 0; (i < curr->ngrp); i++)
for (i = 0; (i < curr->ngrp); i++)
{
/* calculate something useful, like mean square displacements */
- calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != NULL), com,
+ calc_corr(curr, i, gnx[i], index[i], xa[cur], (gnx_com != nullptr), com,
calc1, bTen);
}
cur = prev;
while (read_next_x(oenv, status, &t, x[cur], box));
fprintf(stderr, "\nUsed %d restart points spaced %g %s over %g %s\n\n",
curr->nrestart,
- output_env_conv_time(oenv, dt), output_env_get_time_unit(oenv),
+ output_env_conv_time(oenv, dt), output_env_get_time_unit(oenv).c_str(),
output_env_conv_time(oenv, curr->time[curr->nframes-1]),
- output_env_get_time_unit(oenv) );
+ output_env_get_time_unit(oenv).c_str() );
if (bMol)
{
gmx_rmpbc_done(gpbc);
}
- close_trj(status);
+ close_trx(status);
return natoms;
}
real *DD, *SigmaD, a, a2, b, r, chi2;
rvec *x;
matrix box;
- int *gnx_com = NULL; /* the COM removal group size */
- int **index_com = NULL; /* the COM removal group atom indices */
- char **grpname_com = NULL; /* the COM removal group name */
+ int *gnx_com = nullptr; /* the COM removal group size */
+ int **index_com = nullptr; /* the COM removal group atom indices */
+ char **grpname_com = nullptr; /* the COM removal group name */
snew(gnx, nrgrp);
snew(index, nrgrp);
}
msd = init_corr(nrgrp, type, axis, dim_factor,
- mol_file == NULL ? 0 : gnx[0], bTen, bMW, dt, top,
+ mol_file == nullptr ? 0 : gnx[0], bTen, bMW, dt, top,
beginfit, endfit);
nat_trx =
corr_loop(msd, trx_file, top, ePBC, mol_file ? gnx[0] : 0, gnx, index,
- (mol_file != NULL) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
+ (mol_file != nullptr) ? calc1_mol : (bMW ? calc1_mw : calc1_norm),
bTen, gnx_com, index_com, dt, t_pdb,
- pdb_file ? &x : NULL, box, oenv);
+ pdb_file ? &x : nullptr, box, oenv);
/* Correct for the number of points */
for (j = 0; (j < msd->ngrp); j++)
if (mol_file)
{
- if (pdb_file && x == NULL)
+ if (pdb_file && x == nullptr)
{
fprintf(stderr, "\nNo frame found need time tpdb = %g ps\n"
"Can not write %s\n\n", t_pdb, pdb_file);
}
i = top->atoms.nr;
top->atoms.nr = nat_trx;
- if (pdb_file && top->atoms.pdbinfo == NULL)
+ if (pdb_file && top->atoms.pdbinfo == nullptr)
{
snew(top->atoms.pdbinfo, top->atoms.nr);
}
top->atoms.nr = i;
}
- DD = NULL;
- SigmaD = NULL;
+ DD = nullptr;
+ SigmaD = nullptr;
if (beginfit == -1)
{
}
}
fprintf(stdout, "Fitting from %g to %g %s\n\n", beginfit, endfit,
- output_env_get_time_unit(oenv));
+ output_env_get_time_unit(oenv).c_str());
N = i1-i0;
if (N <= 2)
"The diffusion coefficient is determined by linear regression of the MSD,",
"where, unlike for the normal output of D, the times are weighted",
"according to the number of reference points, i.e. short times have",
- "a higher weight. Also when [TT]-beginfit[tt]=-1,fitting starts at 10%",
- "and when [TT]-endfit[tt]=-1, fitting goes to 90%.",
+ "a higher weight. Also when [TT]-beginfit[tt] is -1, fitting starts at 10%",
+ "and when [TT]-endfit[tt] is -1, fitting goes to 90%.",
"Using this option one also gets an accurate error estimate",
"based on the statistics between individual molecules.",
"Note that this diffusion coefficient and error estimate are only",
"the diffusion coefficient of the molecule.",
"This option implies option [TT]-mol[tt]."
};
- static const char *normtype[] = { NULL, "no", "x", "y", "z", NULL };
- static const char *axtitle[] = { NULL, "no", "x", "y", "z", NULL };
+ static const char *normtype[] = { nullptr, "no", "x", "y", "z", nullptr };
+ static const char *axtitle[] = { nullptr, "no", "x", "y", "z", nullptr };
static int ngroup = 1;
static real dt = 10;
static real t_pdb = 0;
};
t_filenm fnm[] = {
- { efTRX, NULL, NULL, ffREAD },
- { efTPS, NULL, NULL, ffREAD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efXVG, NULL, "msd", ffWRITE },
+ { efTRX, nullptr, nullptr, ffREAD },
+ { efTPS, nullptr, nullptr, ffREAD },
+ { efNDX, nullptr, nullptr, ffOPTRD },
+ { efXVG, nullptr, "msd", ffWRITE },
{ efXVG, "-mol", "diff_mol", ffOPTWR },
{ efPDB, "-pdb", "diff_mol", ffOPTWR }
};
if (!parse_common_args(&argc, argv,
PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_TIME_UNIT,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
fprintf(stderr, "Calculating diffusion coefficients for molecules.\n");
}
- GMX_RELEASE_ASSERT(normtype[0] != 0, "Options inconsistency; normtype[0]Â is NULL");
- GMX_RELEASE_ASSERT(axtitle[0] != 0, "Options inconsistency; axtitle[0]Â is NULL");
+ GMX_RELEASE_ASSERT(normtype[0] != nullptr, "Options inconsistency; normtype[0]Â is NULL");
+ GMX_RELEASE_ASSERT(axtitle[0] != nullptr, "Options inconsistency; axtitle[0]Â is NULL");
if (normtype[0][0] != 'n')
{
gmx_fatal(FARGS, "Can only calculate the full tensor for 3D msd");
}
- bTop = read_tps_conf(tps_file, &top, &ePBC, &xdum, NULL, box, bMW || bRmCOMM);
+ bTop = read_tps_conf(tps_file, &top, &ePBC, &xdum, nullptr, box, bMW || bRmCOMM);
if (mol_file && !bTop)
{
gmx_fatal(FARGS,
#ifndef CUDA_ARCH_UTILS_CUH_
#define CUDA_ARCH_UTILS_CUH_
+ #include "config.h"
+
/*! \file
* \brief CUDA arch dependent definitions.
*
*/
static const int warp_size = 32;
static const int warp_size_log2 = 5;
+ /*! \brief Bitmask corresponding to all threads active in a warp.
+ * NOTE that here too we assume 32-wide warps.
+ */
+ static const unsigned int c_fullWarpMask = 0xffffffff;
+
+ /* Below are backward-compatibility wrappers for CUDA 9 warp-wide intrinsics. */
+
+ /*! \brief Compatibility wrapper around the CUDA __syncwarp() instrinsic. */
+ static __forceinline__ __device__
+ void gmx_syncwarp(const unsigned int activeMask = c_fullWarpMask)
+ {
+ #if GMX_CUDA_VERSION < 9000
+ /* no sync needed on pre-Volta. */
+ GMX_UNUSED_VALUE(activeMask);
+ #else
+ __syncwarp(activeMask);
+ #endif
+ }
+
+ /*! \brief Compatibility wrapper around the CUDA __ballot()/__ballot_sync() instrinsic. */
+ static __forceinline__ __device__
+ unsigned int gmx_ballot_sync(const unsigned int activeMask,
+ const int pred)
+ {
+ #if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __ballot(pred);
+ #else
+ return __ballot_sync(activeMask, pred);
+ #endif
+ }
+
+ /*! \brief Compatibility wrapper around the CUDA __any()/__any_sync() instrinsic. */
+ static __forceinline__ __device__
+ int gmx_any_sync(const unsigned int activeMask,
+ const int pred)
+ {
+ #if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __any(pred);
+ #else
+ return __any_sync(activeMask, pred);
+ #endif
+ }
+
+ /*! \brief Compatibility wrapper around the CUDA __shfl_up()/__shfl_up_sync() instrinsic. */
+ template <typename T>
+ static __forceinline__ __device__
+ T gmx_shfl_up_sync(const unsigned int activeMask,
+ const T var,
+ unsigned int offset)
+ {
+ #if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __shfl_up(var, offset);
+ #else
+ return __shfl_up_sync(activeMask, var, offset);
+ #endif
+ }
+
+ /*! \brief Compatibility wrapper around the CUDA __shfl_down()/__shfl_down_sync() instrinsic. */
+ template <typename T>
+ static __forceinline__ __device__
+ T gmx_shfl_down_sync(const unsigned int activeMask,
+ const T var,
+ unsigned int offset)
+ {
+ #if GMX_CUDA_VERSION < 9000
+ GMX_UNUSED_VALUE(activeMask);
+ return __shfl_down(var, offset);
+ #else
+ return __shfl_down_sync(activeMask, var, offset);
+ #endif
+ }
+/*! \brief Allow disabling CUDA textures using the GMX_DISABLE_CUDA_TEXTURES macro.
+ *
+ * This option will not influence functionality. All features using textures ought
+ * to have fallback for texture-less reads (direct/LDG loads), all new code needs
+ * to provide fallback code.
+ */
+#if defined GMX_DISABLE_CUDA_TEXTURES
+#define DISABLE_CUDA_TEXTURES 1
+#else
+#define DISABLE_CUDA_TEXTURES 0
+#endif
+
+/* CUDA architecture technical characteristics. Needs macros because it is used
+ * in the __launch_bounds__ function qualifiers and might need it in preprocessor
+ * conditionals.
+ *
+ */
+#if GMX_PTX_ARCH > 0
+ #if GMX_PTX_ARCH <= 210 // CC 2.x
+ #define GMX_CUDA_MAX_BLOCKS_PER_MP 8
+ #define GMX_CUDA_MAX_THREADS_PER_MP 1536
+ #elif GMX_PTX_ARCH <= 370 // CC 3.x
+ #define GMX_CUDA_MAX_BLOCKS_PER_MP 16
+ #define GMX_CUDA_MAX_THREADS_PER_MP 2048
+ #else // CC 5.x, 6.x
+ /* Note that this final branch covers all future architectures (current gen
+ * is 6.x as of writing), hence assuming that these *currently defined* upper
+ * limits will not be lowered.
+ */
+ #define GMX_CUDA_MAX_BLOCKS_PER_MP 32
+ #define GMX_CUDA_MAX_THREADS_PER_MP 2048
+ #endif
+#else
+ #define GMX_CUDA_MAX_BLOCKS_PER_MP 0
+ #define GMX_CUDA_MAX_THREADS_PER_MP 0
+#endif
+
#endif /* CUDA_ARCH_UTILS_CUH_ */
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2014,2015,2016, by the GROMACS development team, led by
- * Copyright (c) 2012,2014,2015,2017, by the GROMACS development team, led by
++ * Copyright (c) 2012,2014,2015,2016,2017, 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.
return cu_copy_H2D_generic(d_dest, h_src, bytes, true, s);
}
-int cu_copy_H2D_alloc(void ** d_dest, void * h_src, size_t bytes)
-{
- cudaError_t stat;
-
- if (d_dest == NULL || h_src == NULL || bytes == 0)
- {
- return -1;
- }
-
- stat = cudaMalloc(d_dest, bytes);
- CU_RET_ERR(stat, "cudaMalloc failed in cu_copy_H2D_alloc");
-
- return cu_copy_H2D(*d_dest, h_src, bytes);
-}
-
float cu_event_elapsed(cudaEvent_t start, cudaEvent_t end)
{
float t = 0.0;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2017, 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.
* The CUDA device information is queried and set at detection and contains
* both information about the device/hardware returned by the runtime as well
* as additional data like support status.
+ *
+ * \todo extract an object to manage NVML details
*/
struct gmx_device_info_t
{
int id; /* id of the CUDA device */
cudaDeviceProp prop; /* CUDA device properties */
int stat; /* result of the device check */
- gmx_bool nvml_initialized; /* If NVML was initialized */
unsigned int nvml_orig_app_sm_clock; /* The original SM clock before we changed it */
unsigned int nvml_orig_app_mem_clock; /* The original memory clock before we changed it */
gmx_bool nvml_app_clocks_changed; /* If application clocks have been changed */
unsigned int nvml_set_app_mem_clock; /* The memory clock we set */
#if HAVE_NVML
nvmlDevice_t nvml_device_id; /* NVML device id */
+ // TODO This can become a bool with a more useful name
nvmlEnableState_t nvml_is_restricted; /* Status of application clocks permission */
#endif /* HAVE_NVML */
};
/*! Launches asynchronous host to device memory copy in stream s. */
int cu_copy_D2H_async(void * /*h_dest*/, void * /*d_src*/, size_t /*bytes*/, cudaStream_t /*s = 0*/);
-
/*! Launches synchronous host to device memory copy. */
int cu_copy_H2D(void * /*d_dest*/, void * /*h_src*/, size_t /*bytes*/);
/*! Launches asynchronous host to device memory copy in stream s. */
int cu_copy_H2D_async(void * /*d_dest*/, void * /*h_src*/, size_t /*bytes*/, cudaStream_t /*s = 0*/);
-/*! Allocates device memory and launches synchronous host to device memory copy. */
-int cu_copy_H2D_alloc(void ** /*d_dest*/, void * /*h_src*/, size_t /*bytes*/);
-
/*! Frees device memory and resets the size and allocation size to -1. */
void cu_free_buffered(void *d_ptr, int *n = NULL, int *nalloc = NULL);
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/ewald/ewald.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/mdlib/nbnxn_gpu_data_mgmt.h"
#include "gromacs/mdlib/nbnxn_search.h"
#include "gromacs/mdlib/nbnxn_simd.h"
+#include "gromacs/mdlib/nbnxn_tuning.h"
#include "gromacs/mdlib/nbnxn_util.h"
#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/fcdata.h"
#include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/iforceprovider.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
-#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility/strconvert.h"
#include "nbnxn_gpu_jit_support.h"
const char *egrp_nm[egNR+1] = {
"Coul-SR", "LJ-SR", "Buck-SR",
- "Coul-14", "LJ-14", NULL
+ "Coul-14", "LJ-14", nullptr
};
t_forcerec *mk_forcerec(void)
/* Check if we are doing QM on this group */
qm = FALSE;
- if (qm_grpnr != NULL)
+ if (qm_grpnr != nullptr)
{
for (j = j0; j < j1 && !qm; j++)
{
}
n_solvent_parameters = 0;
- solvent_parameters = NULL;
+ solvent_parameters = nullptr;
/* Allocate temporary array for solvent type */
snew(cg_sp, mtop->nmolblock);
{
check_solvent_cg(molt, cg_mol, nmol,
mtop->groups.grpnr[egcQMMM] ?
- mtop->groups.grpnr[egcQMMM]+at_offset+am : 0,
+ mtop->groups.grpnr[egcQMMM]+at_offset+am : nullptr,
&mtop->groups.grps[egcQMMM],
fr,
&n_solvent_parameters, &solvent_parameters,
}
sfree(cg_sp);
- if (bestsol != esolNO && fp != NULL)
+ if (bestsol != esolNO && fp != nullptr)
{
fprintf(fp, "\nEnabling %s-like water optimization for %d molecules.\n\n",
esol_names[bestsol],
{
bId = FALSE;
}
- if (mtop->groups.grpnr[egcQMMM] != NULL)
+ if (mtop->groups.grpnr[egcQMMM] != nullptr)
{
for (ai = a0; ai < a1; ai++)
{
{
if (fr->eeltype == eelGRF)
{
- calc_rffac(NULL, fr->eeltype, fr->epsilon_r, fr->epsilon_rf,
+ calc_rffac(nullptr, fr->eeltype, fr->epsilon_r, fr->epsilon_rf,
fr->rcoulomb, fr->temp, fr->zsquare, box,
&fr->kappa, &fr->k_rf, &fr->c_rf);
}
int ntp, *typecount;
gmx_bool bBHAM;
real *nbfp;
- real *nbfp_comb = NULL;
+ real *nbfp_comb = nullptr;
ntp = fr->ntype;
bBHAM = fr->bBHAM;
sfree(nbfp_comb);
}
- if (fplog != NULL)
+ if (fplog != nullptr)
{
if (fr->eDispCorr == edispcAllEner ||
fr->eDispCorr == edispcAllEnerPres)
char buf[STRLEN];
int i, j;
- if (tabfn == NULL)
+ if (tabfn == nullptr)
{
if (debug)
{
int ncount, *count;
bondedtable_t *tab;
- tab = NULL;
+ tab = nullptr;
ncount = 0;
- count = NULL;
+ count = nullptr;
count_tables(ftype1, ftype2, mtop, &ncount, &count);
// Are there any relevant tabulated bond interactions?
if (fr->bF_NoVirSum)
{
- fr->f_novirsum_n = natoms_f_novirsum;
- if (fr->f_novirsum_n > fr->f_novirsum_nalloc)
- {
- fr->f_novirsum_nalloc = over_alloc_dd(fr->f_novirsum_n);
- srenew(fr->f_novirsum_alloc, fr->f_novirsum_nalloc);
- }
- }
- else
- {
- fr->f_novirsum_n = 0;
+ /* TODO: remove this + 1 when padding is properly implemented */
+ fr->forceBufferNoVirialSummation->resize(natoms_f_novirsum + 1);
}
}
(ir->implicit_solvent == eisGBSA && (ir->gb_algorithm == egbSTILL ||
ir->gb_algorithm == egbHCT ||
ir->gb_algorithm == egbOBC))) &&
- getenv("GMX_NO_ALLVSALL") == NULL
+ getenv("GMX_NO_ALLVSALL") == nullptr
);
if (bAllvsAll && ir->opts.ngener > 1)
if (bPrintNote)
{
- if (MASTER(cr))
- {
- fprintf(stderr, "\n%s\n", note);
- }
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp, "\n%s\n", note);
}
}
-gmx_bool nbnxn_gpu_acceleration_supported(FILE *fplog,
- const t_commrec *cr,
- const t_inputrec *ir,
- gmx_bool bRerunMD)
-{
- if (bRerunMD && ir->opts.ngener > 1)
- {
- /* Rerun execution time is dominated by I/O and pair search,
- * so GPUs are not very useful, plus they do not support more
- * than one energy group. If the user requested GPUs
- * explicitly, a fatal error is given later. With non-reruns,
- * we fall back to a single whole-of system energy group
- * (which runs much faster than a multiple-energy-groups
- * implementation would), and issue a note in the .log
- * file. Users can re-run if they want the information. */
- md_print_warn(cr, fplog, "Rerun with energy groups is not implemented for GPUs, falling back to the CPU\n");
- return FALSE;
- }
-
- return TRUE;
-}
-
-gmx_bool nbnxn_simd_supported(FILE *fplog,
- const t_commrec *cr,
- const t_inputrec *ir)
+gmx_bool nbnxn_simd_supported(const gmx::MDLogger &mdlog,
+ const t_inputrec *ir)
{
if (ir->vdwtype == evdwPME && ir->ljpme_combination_rule == eljpmeLB)
{
/* LJ PME with LB combination rule does 7 mesh operations.
* This so slow that we don't compile SIMD non-bonded kernels
* for that. */
- md_print_warn(cr, fplog, "LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling back to plain C kernels\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("LJ-PME with Lorentz-Berthelot is not supported with SIMD kernels, falling back to plain C kernels");
return FALSE;
}
#endif /* GMX_NBNXN_SIMD_2XNN && GMX_NBNXN_SIMD_4XN */
- if (getenv("GMX_NBNXN_SIMD_4XN") != NULL)
+ if (getenv("GMX_NBNXN_SIMD_4XN") != nullptr)
{
#ifdef GMX_NBNXN_SIMD_4XN
*kernel_type = nbnxnk4xN_SIMD_4xN;
gmx_fatal(FARGS, "SIMD 4xN kernels requested, but GROMACS has been compiled without support for these kernels");
#endif
}
- if (getenv("GMX_NBNXN_SIMD_2XNN") != NULL)
+ if (getenv("GMX_NBNXN_SIMD_2XNN") != nullptr)
{
#ifdef GMX_NBNXN_SIMD_2XNN
*kernel_type = nbnxnk4xN_SIMD_2xNN;
* With FMA analytical is sometimes faster for a width if 4 as well.
* On BlueGene/Q, this is faster regardless of precision.
* In single precision, this is faster on Bulldozer.
+ * On Skylake table is faster in single and double. TODO: Test 5xxx series.
*/
-#if GMX_SIMD_REAL_WIDTH >= 8 || \
- (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE) || GMX_SIMD_IBM_QPX
+#if ((GMX_SIMD_REAL_WIDTH >= 8 || (GMX_SIMD_REAL_WIDTH >= 4 && GMX_SIMD_HAVE_FMA && !GMX_DOUBLE)) \
+ && !GMX_SIMD_X86_AVX_512) || GMX_SIMD_IBM_QPX
*ewald_excl = ewaldexclAnalytical;
#endif
- if (getenv("GMX_NBNXN_EWALD_TABLE") != NULL)
+ if (getenv("GMX_NBNXN_EWALD_TABLE") != nullptr)
{
*ewald_excl = ewaldexclTable;
}
- if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != NULL)
+ if (getenv("GMX_NBNXN_EWALD_ANALYTICAL") != nullptr)
{
*ewald_excl = ewaldexclAnalytical;
}
const char *lookup_nbnxn_kernel_name(int kernel_type)
{
- const char *returnvalue = NULL;
+ const char *returnvalue = nullptr;
switch (kernel_type)
{
case nbnxnkNotSet:
case nbnxnkNR:
default:
gmx_fatal(FARGS, "Illegal kernel type selected");
- returnvalue = NULL;
+ returnvalue = nullptr;
break;
}
return returnvalue;
};
static void pick_nbnxn_kernel(FILE *fp,
- const t_commrec *cr,
+ const gmx::MDLogger &mdlog,
gmx_bool use_simd_kernels,
gmx_bool bUseGPU,
- gmx_bool bEmulateGPU,
+ bool emulateGpu,
const t_inputrec *ir,
int *kernel_type,
int *ewald_excl,
*kernel_type = nbnxnkNotSet;
*ewald_excl = ewaldexclTable;
- if (bEmulateGPU)
+ if (emulateGpu)
{
*kernel_type = nbnxnk8x8x8_PlainC;
if (bDoNonbonded)
{
- md_print_warn(cr, fp, "Emulating a GPU run on the CPU (slow)");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("Emulating a GPU run on the CPU (slow)");
}
}
else if (bUseGPU)
if (*kernel_type == nbnxnkNotSet)
{
if (use_simd_kernels &&
- nbnxn_simd_supported(fp, cr, ir))
+ nbnxn_simd_supported(mdlog, ir))
{
pick_nbnxn_kernel_cpu(ir, kernel_type, ewald_excl);
}
}
}
- if (bDoNonbonded && fp != NULL)
+ if (bDoNonbonded && fp != nullptr)
{
fprintf(fp, "\nUsing %s %dx%d non-bonded kernels\n\n",
lookup_nbnxn_kernel_name(*kernel_type),
if (nbnxnk4x4_PlainC == *kernel_type ||
nbnxnk8x8x8_PlainC == *kernel_type)
{
- md_print_warn(cr, fp,
- "WARNING: Using the slow %s kernels. This should\n"
- "not happen during routine usage on supported platforms.\n\n",
- lookup_nbnxn_kernel_name(*kernel_type));
+ GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+ "WARNING: Using the slow %s kernels. This should\n"
+ "not happen during routine usage on supported platforms.",
+ lookup_nbnxn_kernel_name(*kernel_type));
}
}
}
-static void pick_nbnxn_resources(FILE *fp,
- const t_commrec *cr,
- const gmx_hw_info_t *hwinfo,
- gmx_bool bDoNonbonded,
- gmx_bool *bUseGPU,
- gmx_bool *bEmulateGPU,
- const gmx_gpu_opt_t *gpu_opt)
-{
- gmx_bool bEmulateGPUEnvVarSet;
- char gpu_err_str[STRLEN];
-
- *bUseGPU = FALSE;
-
- bEmulateGPUEnvVarSet = (getenv("GMX_EMULATE_GPU") != NULL);
-
- /* Run GPU emulation mode if GMX_EMULATE_GPU is defined. Because
- * GPUs (currently) only handle non-bonded calculations, we will
- * automatically switch to emulation if non-bonded calculations are
- * turned off via GMX_NO_NONBONDED - this is the simple and elegant
- * way to turn off GPU initialization, data movement, and cleanup.
- *
- * GPU emulation can be useful to assess the performance one can expect by
- * adding GPU(s) to the machine. The conditional below allows this even
- * if mdrun is compiled without GPU acceleration support.
- * Note that you should freezing the system as otherwise it will explode.
- */
- *bEmulateGPU = (bEmulateGPUEnvVarSet ||
- (!bDoNonbonded && gpu_opt->n_dev_use > 0));
-
- /* Enable GPU mode when GPUs are available or no GPU emulation is requested.
- */
- if (gpu_opt->n_dev_use > 0 && !(*bEmulateGPU))
- {
- /* Each PP node will use the intra-node id-th device from the
- * list of detected/selected GPUs. */
- if (!init_gpu(fp, cr->rank_pp_intranode, gpu_err_str,
- &hwinfo->gpu_info, gpu_opt))
- {
- /* At this point the init should never fail as we made sure that
- * we have all the GPUs we need. If it still does, we'll bail. */
- /* TODO the decorating of gpu_err_str is nicer if it
- happens inside init_gpu. Out here, the decorating with
- the MPI rank makes sense. */
- gmx_fatal(FARGS, "On rank %d failed to initialize GPU #%d: %s",
- cr->nodeid,
- get_gpu_device_id(&hwinfo->gpu_info, gpu_opt,
- cr->rank_pp_intranode),
- gpu_err_str);
- }
-
- /* Here we actually turn on hardware GPU acceleration */
- *bUseGPU = TRUE;
- }
-}
-
gmx_bool uses_simple_tables(int cutoff_scheme,
nonbonded_verlet_t *nbv,
int group)
{
init_ewald_f_table(ic, rtab);
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp, "Initialized non-bonded Ewald correction tables, spacing: %.2e size: %d\n\n",
1/ic->tabq_scale, ic->tabq_size);
snew_aligned(ic->tabq_coul_F, 16, 32);
snew_aligned(ic->tabq_coul_V, 16, 32);
- ic->rlist = fr->rlist;
-
/* Lennard-Jones */
ic->vdwtype = fr->vdwtype;
ic->vdw_modifier = fr->vdw_modifier;
ic->epsfac = fr->epsfac;
ic->ewaldcoeff_q = fr->ewaldcoeff_q;
- if (fr->coulomb_modifier == eintmodPOTSHIFT)
+ if (EEL_PME_EWALD(ic->eeltype) && ic->coulomb_modifier == eintmodPOTSHIFT)
{
- ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
+ GMX_RELEASE_ASSERT(ic->rcoulomb != 0, "Cutoff radius cannot be zero");
+ ic->sh_ewald = std::erfc(ic->ewaldcoeff_q*ic->rcoulomb) / ic->rcoulomb;
}
else
{
}
}
- if (fp != NULL)
+ if (fp != nullptr)
{
real dispersion_shift;
*interaction_const = ic;
}
+/* TODO deviceInfo should be logically const, but currently
+ * init_gpu modifies it to set up NVML support. This could
+ * happen during the detection phase, and deviceInfo could
+ * the become const. */
static void init_nb_verlet(FILE *fp,
+ const gmx::MDLogger &mdlog,
nonbonded_verlet_t **nb_verlet,
gmx_bool bFEP_NonBonded,
const t_inputrec *ir,
const t_forcerec *fr,
const t_commrec *cr,
- const char *nbpu_opt)
+ const char *nbpu_opt,
+ gmx_device_info_t *deviceInfo,
+ const gmx_mtop_t *mtop,
+ matrix box)
{
nonbonded_verlet_t *nbv;
int i;
char *env;
- gmx_bool bEmulateGPU, bHybridGPURun = FALSE;
+ gmx_bool bHybridGPURun = FALSE;
nbnxn_alloc_t *nb_alloc;
nbnxn_free_t *nb_free;
- snew(nbv, 1);
+ nbv = new nonbonded_verlet_t();
- pick_nbnxn_resources(fp, cr, fr->hwinfo,
- fr->bNonbonded,
- &nbv->bUseGPU,
- &bEmulateGPU,
- fr->gpu_opt);
+ nbv->emulateGpu = (getenv("GMX_EMULATE_GPU") != nullptr);
+ nbv->bUseGPU = deviceInfo != nullptr;
- nbv->nbs = NULL;
+ GMX_RELEASE_ASSERT(!(nbv->emulateGpu && nbv->bUseGPU), "When GPU emulation is active, there cannot be a GPU assignment");
+
+ if (nbv->bUseGPU)
+ {
+ /* Use the assigned GPU. */
+ init_gpu(mdlog, cr->nodeid, deviceInfo);
+ }
+
+ nbv->nbs = nullptr;
nbv->min_ci_balanced = 0;
nbv->ngrp = (DOMAINDECOMP(cr) ? 2 : 1);
for (i = 0; i < nbv->ngrp; i++)
{
nbv->grp[i].nbl_lists.nnbl = 0;
- nbv->grp[i].nbat = NULL;
+ nbv->grp[i].nbat = nullptr;
nbv->grp[i].kernel_type = nbnxnkNotSet;
if (i == 0) /* local */
{
- pick_nbnxn_kernel(fp, cr, fr->use_simd_kernels,
- nbv->bUseGPU, bEmulateGPU, ir,
+ pick_nbnxn_kernel(fp, mdlog, fr->use_simd_kernels,
+ nbv->bUseGPU, nbv->emulateGpu, ir,
&nbv->grp[i].kernel_type,
&nbv->grp[i].ewald_excl,
fr->bNonbonded);
}
else /* non-local */
{
- if (nbpu_opt != NULL && strcmp(nbpu_opt, "gpu_cpu") == 0)
+ if (nbpu_opt != nullptr && strcmp(nbpu_opt, "gpu_cpu") == 0)
{
/* Use GPU for local, select a CPU kernel for non-local */
- pick_nbnxn_kernel(fp, cr, fr->use_simd_kernels,
- FALSE, FALSE, ir,
+ pick_nbnxn_kernel(fp, mdlog, fr->use_simd_kernels,
+ FALSE, false, ir,
&nbv->grp[i].kernel_type,
&nbv->grp[i].ewald_excl,
fr->bNonbonded);
}
}
+ nbv->listParams = std::unique_ptr<NbnxnListParameters>(new NbnxnListParameters(ir->rlist));
+ setupDynamicPairlistPruning(fp, ir, mtop, box, nbv->bUseGPU, fr->ic,
+ nbv->listParams.get());
+
nbnxn_init_search(&nbv->nbs,
- DOMAINDECOMP(cr) ? &cr->dd->nc : NULL,
- DOMAINDECOMP(cr) ? domdec_zones(cr->dd) : NULL,
+ DOMAINDECOMP(cr) ? &cr->dd->nc : nullptr,
+ DOMAINDECOMP(cr) ? domdec_zones(cr->dd) : nullptr,
bFEP_NonBonded,
gmx_omp_nthreads_get(emntPairsearch));
if (fr->vdwtype == evdwCUT &&
(fr->vdw_modifier == eintmodNONE ||
fr->vdw_modifier == eintmodPOTSHIFT) &&
- getenv("GMX_NO_LJ_COMB_RULE") == NULL)
+ getenv("GMX_NO_LJ_COMB_RULE") == nullptr)
{
/* Plain LJ cut-off: we can optimize with combination rules */
enbnxninitcombrule = enbnxninitcombruleDETECT;
/* init the NxN GPU data; the last argument tells whether we'll have
* both local and non-local NB calculation on GPU */
nbnxn_gpu_init(&nbv->gpu_nbv,
- &fr->hwinfo->gpu_info,
- fr->gpu_opt,
+ deviceInfo,
fr->ic,
+ nbv->listParams.get(),
nbv->grp,
- cr->rank_pp_intranode,
cr->nodeid,
(nbv->ngrp > 1) && !bHybridGPURun);
}
#endif /* GMX_THREAD_MPI */
- if ((env = getenv("GMX_NB_MIN_CI")) != NULL)
+ if ((env = getenv("GMX_NB_MIN_CI")) != nullptr)
{
char *end;
gmx_bool usingGpu(nonbonded_verlet_t *nbv)
{
- return nbv != NULL && nbv->bUseGPU;
+ return nbv != nullptr && nbv->bUseGPU;
}
-void init_forcerec(FILE *fp,
- t_forcerec *fr,
- t_fcdata *fcd,
- const t_inputrec *ir,
- const gmx_mtop_t *mtop,
- const t_commrec *cr,
- matrix box,
- const char *tabfn,
- const char *tabpfn,
- const t_filenm *tabbfnm,
- const char *nbpu_opt,
- gmx_bool bNoSolvOpt,
- real print_force)
+void init_forcerec(FILE *fp,
+ const gmx::MDLogger &mdlog,
+ t_forcerec *fr,
+ t_fcdata *fcd,
+ const t_inputrec *ir,
+ const gmx_mtop_t *mtop,
+ const t_commrec *cr,
+ matrix box,
+ const char *tabfn,
+ const char *tabpfn,
+ const t_filenm *tabbfnm,
+ const char *nbpu_opt,
+ gmx_device_info_t *deviceInfo,
+ gmx_bool bNoSolvOpt,
+ real print_force)
{
int i, m, negp_pp, negptable, egi, egj;
real rtab;
gmx_bool bFEP_NonBonded;
int *nm_ind, egp_flags;
- if (fr->hwinfo == NULL)
- {
- /* Detect hardware, gather information.
- * In mdrun, hwinfo has already been set before calling init_forcerec.
- * Here we ignore GPUs, as tools will not use them anyhow.
- */
- fr->hwinfo = gmx_detect_hardware(fp, cr, FALSE);
- }
-
/* By default we turn SIMD kernels on, but it might be turned off further down... */
fr->use_simd_kernels = TRUE;
fr->sc_sigma6_def = gmx::power6(ir->fepvals->sc_sigma);
env = getenv("GMX_SCSIGMA_MIN");
- if (env != NULL)
+ if (env != nullptr)
{
dbl = 0;
sscanf(env, "%20lf", &dbl);
}
fr->bNonbonded = TRUE;
- if (getenv("GMX_NO_NONBONDED") != NULL)
+ if (getenv("GMX_NO_NONBONDED") != nullptr)
{
/* turn off non-bonded calculations */
fr->bNonbonded = FALSE;
- md_print_warn(cr, fp,
- "Found environment variable GMX_NO_NONBONDED.\n"
- "Disabling nonbonded calculations.\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText(
+ "Found environment variable GMX_NO_NONBONDED.\n"
+ "Disabling nonbonded calculations.");
}
bGenericKernelOnly = FALSE;
* can be used with water optimization, and disable it if that is not the case.
*/
- if (getenv("GMX_NB_GENERIC") != NULL)
+ if (getenv("GMX_NB_GENERIC") != nullptr)
{
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp,
"Found environment variable GMX_NB_GENERIC.\n"
bNoSolvOpt = TRUE;
}
- if ( (getenv("GMX_DISABLE_SIMD_KERNELS") != NULL) || (getenv("GMX_NOOPTIMIZEDKERNELS") != NULL) )
+ if ( (getenv("GMX_DISABLE_SIMD_KERNELS") != nullptr) || (getenv("GMX_NOOPTIMIZEDKERNELS") != nullptr) )
{
fr->use_simd_kernels = FALSE;
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp,
"\nFound environment variable GMX_DISABLE_SIMD_KERNELS.\n"
fr->bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
/* Check if we can/should do all-vs-all kernels */
- fr->bAllvsAll = can_use_allvsall(ir, FALSE, NULL, NULL);
- fr->AllvsAll_work = NULL;
- fr->AllvsAll_workgb = NULL;
+ fr->bAllvsAll = can_use_allvsall(ir, FALSE, nullptr, nullptr);
+ fr->AllvsAll_work = nullptr;
+ fr->AllvsAll_workgb = nullptr;
/* All-vs-all kernels have not been implemented in 4.6 and later.
* See Redmine #1249. */
if (fr->bAllvsAll)
{
fr->bAllvsAll = FALSE;
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp,
"\nYour simulation settings would have triggered the efficient all-vs-all\n"
{
fprintf(stderr, "\n%s\n", note);
}
- if (fp != NULL)
+ if (fp != nullptr)
{
fprintf(fp, "\n%s\n", note);
}
}
else
{
- fr->bMolPBC = TRUE;
+ /* Not making molecules whole is faster in most cases,
+ * but With orientation restraints we need whole molecules.
+ */
+ fr->bMolPBC = (fcd->orires.nr == 0);
- if (getenv("GMX_USE_GRAPH") != NULL)
+ if (getenv("GMX_USE_GRAPH") != nullptr)
{
fr->bMolPBC = FALSE;
if (fp)
{
- md_print_warn(cr, fp, "GMX_USE_GRAPH is set, using the graph for bonded interactions\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("GMX_USE_GRAPH is set, using the graph for bonded interactions");
}
if (mtop->bIntermolecularInteractions)
{
- md_print_warn(cr, fp, "WARNING: Molecules linked by intermolecular interactions have to reside in the same periodic image, otherwise artifacts will occur!\n");
+ GMX_LOG(mdlog.warning).asParagraph().appendText("WARNING: Molecules linked by intermolecular interactions have to reside in the same periodic image, otherwise artifacts will occur!");
}
}
+ GMX_RELEASE_ASSERT(fr->bMolPBC || !mtop->bIntermolecularInteractions, "We need to use PBC within molecules with inter-molecular interactions");
+
if (bSHAKE && fr->bMolPBC)
{
gmx_fatal(FARGS, "SHAKE is not properly supported with intermolecular interactions. For short simulations where linked molecules remain in the same periodic image, the environment variable GMX_USE_GRAPH can be used to override this check.\n");
fr->bcoultab = FALSE;
}
+ /* This now calculates sum for q and C6 */
+ set_chargesum(fp, fr, mtop);
+
/* Tables are used for direct ewald sum */
if (fr->bEwald)
{
if (ir->ewald_geometry == eewg3DC)
{
+ bool haveNetCharge = (fabs(fr->qsum[0]) > 1e-4 ||
+ fabs(fr->qsum[1]) > 1e-4);
if (fp)
{
- fprintf(fp, "Using the Ewald3DC correction for systems with a slab geometry.\n");
+ fprintf(fp, "Using the Ewald3DC correction for systems with a slab geometry%s.\n",
+ haveNetCharge ? " and net charge" : "");
}
please_cite(fp, "In-Chul99a");
+ if (haveNetCharge)
+ {
+ please_cite(fp, "Ballenegger2009");
+ }
}
}
fr->ewaldcoeff_q = calc_ewaldcoeff_q(ir->rcoulomb, ir->ewald_rtol);
}
fr->bF_NoVirSum = (EEL_FULL(fr->eeltype) || EVDW_PME(fr->vdwtype) ||
+ fr->forceProviders->hasForcesWithoutVirialContribution() ||
gmx_mtop_ftype_count(mtop, F_POSRES) > 0 ||
- gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0 ||
- inputrecElecField(ir)
- );
+ gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0);
+
+ if (fr->bF_NoVirSum)
+ {
+ fr->forceBufferNoVirialSummation = new PaddedRVecVector;
+ }
if (fr->cutoff_scheme == ecutsGROUP &&
ncg_mtop(mtop) > fr->cg_nalloc && !DOMAINDECOMP(cr))
fr->cg_nalloc = ncg_mtop(mtop);
srenew(fr->cg_cm, fr->cg_nalloc);
}
- if (fr->shift_vec == NULL)
+ if (fr->shift_vec == nullptr)
{
snew(fr->shift_vec, SHIFTS);
}
- if (fr->fshift == NULL)
+ if (fr->fshift == nullptr)
{
snew(fr->fshift, SHIFTS);
}
- if (fr->nbfp == NULL)
+ if (fr->nbfp == nullptr)
{
fr->ntype = mtop->ffparams.atnr;
fr->nbfp = mk_nbfp(&mtop->ffparams, fr->bBHAM);
&fr->kappa, &fr->k_rf, &fr->c_rf);
}
- /*This now calculates sum for q and c6*/
- set_chargesum(fp, fr, mtop);
-
/* Construct tables for the group scheme. A little unnecessary to
* make both vdw and coul tables sometimes, but what the
* heck. Note that both cutoff schemes construct Ewald tables in
/* make tables for ordinary interactions */
if (bSomeNormalNbListsAreInUse)
{
- make_nbf_tables(fp, fr, rtab, tabfn, NULL, NULL, &fr->nblists[0]);
+ make_nbf_tables(fp, fr, rtab, tabfn, nullptr, nullptr, &fr->nblists[0]);
m = 1;
}
else
&fr->bExcl_IntraCGAll_InterCGNone);
if (DOMAINDECOMP(cr))
{
- fr->cginfo = NULL;
+ fr->cginfo = nullptr;
}
else
{
GMX_RELEASE_ASSERT(ir->rcoulomb == ir->rvdw, "With Verlet lists and no PME rcoulomb and rvdw should be identical");
}
- init_nb_verlet(fp, &fr->nbv, bFEP_NonBonded, ir, fr, cr, nbpu_opt);
+ init_nb_verlet(fp, mdlog, &fr->nbv, bFEP_NonBonded, ir, fr,
+ cr, nbpu_opt, deviceInfo,
+ mtop, box);
}
if (ir->eDispCorr != edispcNO)
* in this run because the PME ranks have no knowledge of whether GPUs
* are used or not, but all ranks need to enter the barrier below.
*/
-void free_gpu_resources(const t_forcerec *fr,
- const t_commrec *cr,
- const gmx_gpu_info_t *gpu_info,
- const gmx_gpu_opt_t *gpu_opt)
+void free_gpu_resources(const t_forcerec *fr,
+ const t_commrec *cr,
+ const gmx_device_info_t *deviceInfo)
{
gmx_bool bIsPPrankUsingGPU;
char gpu_err_str[STRLEN];
nbnxn_gpu_free(fr->nbv->gpu_nbv);
/* stop the GPU profiler (only CUDA) */
stopGpuProfiler();
+ }
- /* With tMPI we need to wait for all ranks to finish deallocation before
- * destroying the CUDA context in free_gpu() as some tMPI ranks may be sharing
- * GPU and context.
- *
- * This is not a concern in OpenCL where we use one context per rank which
- * is freed in nbnxn_gpu_free().
- *
- * Note: as only PP ranks need to free GPU resources, so it is safe to
- * not call the barrier on PME ranks.
- */
+ /* With tMPI we need to wait for all ranks to finish deallocation before
+ * destroying the CUDA context in free_gpu() as some tMPI ranks may be sharing
+ * GPU and context.
+ *
+ * This is not a concern in OpenCL where we use one context per rank which
+ * is freed in nbnxn_gpu_free().
+ *
+ * Note: it is safe to not call the barrier on the ranks which do not use GPU,
+ * but it is easier and more futureproof to call it on the whole node.
+ */
#if GMX_THREAD_MPI
- if (PAR(cr))
- {
- gmx_barrier(cr);
- }
+ if (PAR(cr) || MULTISIM(cr))
+ {
+ gmx_barrier_physical_node(cr);
+ }
#endif /* GMX_THREAD_MPI */
+ if (bIsPPrankUsingGPU)
+ {
/* uninitialize GPU (by destroying the context) */
- if (!free_cuda_gpu(cr->rank_pp_intranode, gpu_err_str, gpu_info, gpu_opt))
+ if (!free_cuda_gpu(deviceInfo, gpu_err_str))
{
gmx_warning("On rank %d failed to free GPU #%d: %s",
cr->nodeid, get_current_cuda_gpu_device_id(), gpu_err_str);
*/
#include "gromacs/gpu_utils/cuda_arch_utils.cuh"
+#include "gromacs/gpu_utils/cuda_kernel_utils.cuh"
#include "gromacs/math/utilities.h"
#include "gromacs/pbcutil/ishift.h"
/* Note that floating-point constants in CUDA code should be suffixed
* code that is in double precision.
*/
-#if GMX_PTX_ARCH < 300
-#error "nbnxn_cuda_kernel.cuh included with GMX_PTX_ARCH < 300"
+#if GMX_PTX_ARCH < 300 && GMX_PTX_ARCH != 0
+#error "nbnxn_cuda_kernel.cuh included with GMX_PTX_ARCH < 300 or host pass"
#endif
#if defined EL_EWALD_ANA || defined EL_EWALD_TAB
Kernel launch parameters:
- #blocks = #pair lists, blockId = pair list Id
- #threads = NTHREAD_Z * c_clSize^2
- - shmem = see nbnxn_cuda.cu:calc_shmem_required()
+ - shmem = see nbnxn_cuda.cu:calc_shmem_required_nonbonded()
Each thread calculates an i force-component taking one pair of i-j atoms.
*/
*
* Note that the current kernel implementation only supports NTHREAD_Z > 1 with
* shuffle-based reduction, hence CC >= 3.0.
+ *
+ *
+ * NOTEs / TODO on Volta / CUDA 9 support extensions:
+ * - the current way of computing active mask using ballot_sync() should be
+ * reconsidered: we can compute all masks with bitwise ops iso ballot and
+ * secondly, all conditionals are warp-uniform, so the sync is not needed;
+ * - reconsider the use of __syncwarp(): its only role is currently to prevent
+ * WAR hazard due to the cj preload; we should try to replace it with direct
+ * loads (which may be faster given the improved L1 on Volta).
*/
/* Kernel launch bounds for different compute capabilities. The value of NTHREAD_Z
#ifdef EL_RF
float two_k_rf = nbparam.two_k_rf;
#endif
-#ifdef EL_EWALD_TAB
- float coulomb_tab_scale = nbparam.coulomb_tab_scale;
-#endif
#ifdef EL_EWALD_ANA
float beta2 = nbparam.ewald_beta*nbparam.ewald_beta;
float beta3 = nbparam.ewald_beta*nbparam.ewald_beta*nbparam.ewald_beta;
#endif
#ifdef PRUNE_NBL
- float rlist_sq = nbparam.rlist_sq;
+ float rlist_sq = nbparam.rlistOuter_sq;
#endif
#ifdef CALC_ENERGIES
ci = sci * c_numClPerSupercl + tidxj;
ai = ci * c_clSize + tidxi;
- xqbuf = xq[ai] + shift_vec[nb_sci.shift];
+ float *shiftptr = (float *)&shift_vec[nb_sci.shift];
+ xqbuf = xq[ai] + make_float4(LDG(shiftptr), LDG(shiftptr + 1), LDG(shiftptr + 2), 0.0f);
xqbuf.w *= nbparam.epsfac;
xqib[tidxj * c_clSize + tidxi] = xqbuf;
#endif
#ifdef LJ_EWALD
+ #if DISABLE_CUDA_TEXTURES
+ E_lj += LDG(&nbparam.nbfp[atom_types[(sci*c_numClPerSupercl + i)*c_clSize + tidxi]*(ntypes + 1)*2]);
+ #else
E_lj += tex1Dfetch<float>(nbparam.nbfp_texobj, atom_types[(sci*c_numClPerSupercl + i)*c_clSize + tidxi]*(ntypes + 1)*2);
+ #endif
#endif
}
#endif /* CALC_ENERGIES */
- const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
+#ifdef EXCLUSION_FORCES
++ const int nonSelfInteraction = !(nb_sci.shift == CENTRAL & tidxj <= tidxi);
+#endif
+
+ int j4LoopStart = cij4_start + tidxz;
+ unsigned int j4LoopThreadMask = gmx_ballot_sync(c_fullWarpMask, j4LoopStart < cij4_end);
/* loop over the j clusters = seen by any of the atoms in the current super-cluster */
- for (j4 = cij4_start + tidxz; j4 < cij4_end; j4 += NTHREAD_Z)
+ for (j4 = j4LoopStart; j4 < cij4_end; j4 += NTHREAD_Z)
{
wexcl_idx = pl_cj4[j4].imei[widx].excl_ind;
imask = pl_cj4[j4].imei[widx].imask;
wexcl = excl[wexcl_idx].pair[(tidx) & (warp_size - 1)];
+ unsigned int imaskSkipConditionThreadMask = j4LoopThreadMask;
#ifndef PRUNE_NBL
+ imaskSkipConditionThreadMask = gmx_ballot_sync(j4LoopThreadMask, imask);
if (imask)
#endif
{
/* Pre-load cj into shared memory on both warps separately */
- if ((tidxj == 0 || tidxj == 4) && tidxi < c_nbnxnGpuJgroupSize)
+ if ((tidxj == 0 | tidxj == 4) & (tidxi < c_nbnxnGpuJgroupSize))
{
cjs[tidxi + tidxj * c_nbnxnGpuJgroupSize/c_splitClSize] = pl_cj4[j4].cj[tidxi];
}
+ gmx_syncwarp(imaskSkipConditionThreadMask);
/* Unrolling this loop
- with pruning leads to register spilling;
Tested with up to nvcc 7.5 */
for (jm = 0; jm < c_nbnxnGpuJgroupSize; jm++)
{
- if (imask & (superClInteractionMask << (jm * c_numClPerSupercl)))
+ const unsigned int jmSkipCondition = imask & (superClInteractionMask << (jm * c_numClPerSupercl));
+ const unsigned int jmSkipConditionThreadMask = gmx_ballot_sync(imaskSkipConditionThreadMask, jmSkipCondition);
+ if (jmSkipCondition)
{
mask_ji = (1U << (jm * c_numClPerSupercl));
#endif
for (i = 0; i < c_numClPerSupercl; i++)
{
- if (imask & mask_ji)
+ const unsigned int iInnerSkipCondition = imask & mask_ji;
+ const unsigned int iInnerSkipConditionThreadMask = gmx_ballot_sync(jmSkipConditionThreadMask, iInnerSkipCondition);
+ if (iInnerSkipCondition)
{
ci = sci * c_numClPerSupercl + i; /* i cluster index */
/* If _none_ of the atoms pairs are in cutoff range,
the bit corresponding to the current
cluster-pair in imask gets set to 0. */
- if (!__any(r2 < rlist_sq))
+ if (!gmx_any_sync(iInnerSkipConditionThreadMask, r2 < rlist_sq))
{
imask &= ~mask_ji;
}
/* cutoff & exclusion check */
#ifdef EXCLUSION_FORCES
- if (r2 < rcoulomb_sq *
- (nb_sci.shift != CENTRAL || ci != cj || tidxj > tidxi))
+ if ((r2 < rcoulomb_sq) * (nonSelfInteraction | (ci != cj)))
#else
- if (r2 < rcoulomb_sq * int_bit)
+ if ((r2 < rcoulomb_sq) * int_bit)
#endif
{
/* load the rest of the i-atom parameters */
#ifndef LJ_COMB
/* LJ 6*C6 and 12*C12 */
typei = atib[i * c_clSize + tidxi];
- c6 = tex1Dfetch<float>(nbparam.nbfp_texobj, 2 * (ntypes * typei + typej));
- c12 = tex1Dfetch<float>(nbparam.nbfp_texobj, 2 * (ntypes * typei + typej) + 1);
+ fetch_nbfp_c6_c12(c6, c12, nbparam, ntypes * typei + typej);
#else
ljcp_i = ljcpib[i * c_clSize + tidxi];
#ifdef LJ_COMB_GEOM
#ifdef LJ_POT_SWITCH
#ifdef CALC_ENERGIES
- calculate_potential_switch_F_E(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+ calculate_potential_switch_F_E(nbparam, inv_r, r2, &F_invr, &E_lj_p);
#else
- calculate_potential_switch_F(nbparam, c6, c12, inv_r, r2, &F_invr, &E_lj_p);
+ calculate_potential_switch_F(nbparam, inv_r, r2, &F_invr, &E_lj_p);
#endif /* CALC_ENERGIES */
#endif /* LJ_POT_SWITCH */
F_invr += qi * qj_f * (int_bit*inv_r2*inv_r + pmecorrF(beta2*r2)*beta3);
#elif defined EL_EWALD_TAB
F_invr += qi * qj_f * (int_bit*inv_r2 -
- interpolate_coulomb_force_r(nbparam.coulomb_tab_texobj, r2 * inv_r, coulomb_tab_scale)) * inv_r;
+ interpolate_coulomb_force_r(nbparam, r2 * inv_r)) * inv_r;
#endif /* EL_EWALD_ANA/TAB */
#ifdef CALC_ENERGIES
}
/* reduce j forces */
- reduce_force_j_warp_shfl(fcj_buf, f, tidxi, aj);
+ reduce_force_j_warp_shfl(fcj_buf, f, tidxi, aj, jmSkipConditionThreadMask);
}
}
#ifdef PRUNE_NBL
pl_cj4[j4].imei[widx].imask = imask;
#endif
}
+ // avoid shared memory WAR hazards between loop iterations
+ gmx_syncwarp(j4LoopThreadMask);
+ // update thread mask for next loop iteration
+ j4LoopThreadMask = gmx_ballot_sync(j4LoopThreadMask, (j4 + NTHREAD_Z) < cij4_end);
}
/* skip central shifts when summing shift forces */
ai = (sci * c_numClPerSupercl + i) * c_clSize + tidxi;
reduce_force_i_warp_shfl(fci_buf[i], f,
&fshift_buf, bCalcFshift,
- tidxj, ai);
+ tidxj, ai, c_fullWarpMask);
}
/* add up local shift forces into global mem, tidxj indexes x,y,z */
#ifdef CALC_ENERGIES
/* reduce the energies over warps and store into global memory */
- reduce_energy_warp_shfl(E_lj, E_el, e_lj, e_el, tidx);
+ reduce_energy_warp_shfl(E_lj, E_el, e_lj, e_el, tidx, c_fullWarpMask);
#endif
}
#endif /* FUNCTION_DECLARATION_ONLY */
* \author Szilárd Páll <pall.szilard@gmail.com>
* \ingroup module_mdlib
*/
-#include "config.h"
-
#include <assert.h>
/* Note that floating-point constants in CUDA code should be suffixed
* with f (e.g. 0.5f), to stop the compiler producing intermediate
* code that is in double precision.
*/
-#include "config.h"
#include "gromacs/gpu_utils/cuda_arch_utils.cuh"
+#include "gromacs/gpu_utils/cuda_kernel_utils.cuh"
#include "gromacs/gpu_utils/vectype_ops.cuh"
#include "nbnxn_cuda_types.h"
#ifndef NBNXN_CUDA_KERNEL_UTILS_CUH
#define NBNXN_CUDA_KERNEL_UTILS_CUH
-/* Use texture objects if supported by the target hardware. */
-#if GMX_PTX_ARCH >= 300
+/* Use texture objects if supported by the target hardware (and in host pass). */
+#if GMX_PTX_ARCH >= 300 || GMX_PTX_ARCH == 0
/* Note: convenience macro, needs to be undef-ed at the end of the file. */
#define USE_TEXOBJ
#endif
/*! \brief Log of the i and j cluster size.
* change this together with c_clSize !*/
- static const int c_clSizeLog2 = 3;
+ static const int c_clSizeLog2 = 3;
/*! \brief Square of cluster size. */
- static const int c_clSizeSq = c_clSize*c_clSize;
+ static const int c_clSizeSq = c_clSize*c_clSize;
/*! \brief j-cluster size after split (4 in the current implementation). */
- static const int c_splitClSize = c_clSize/c_nbnxnGpuClusterpairSplit;
+ static const int c_splitClSize = c_clSize/c_nbnxnGpuClusterpairSplit;
/*! \brief Stride in the force accumualation buffer */
- static const int c_fbufStride = c_clSizeSq;
+ static const int c_fbufStride = c_clSizeSq;
+/*! \brief i-cluster interaction mask for a super-cluster with all c_numClPerSupercl=8 bits set */
- static const unsigned superClInteractionMask = ((1U << c_numClPerSupercl) - 1U);
++static const unsigned superClInteractionMask = ((1U << c_numClPerSupercl) - 1U);
- static const float c_oneSixth = 0.16666667f;
- static const float c_oneTwelveth = 0.08333333f;
+ static const float c_oneSixth = 0.16666667f;
+ static const float c_oneTwelveth = 0.08333333f;
-/* With multiple compilation units this ensures that texture refs are available
- in the the kernels' compilation units. */
-#if !GMX_CUDA_NB_SINGLE_COMPILATION_UNIT
-/*! Texture reference for LJ C6/C12 parameters; bound to cu_nbparam_t.nbfp */
-extern texture<float, 1, cudaReadModeElementType> nbfp_texref;
-
-/*! Texture reference for LJ-PME parameters; bound to cu_nbparam_t.nbfp_comb */
-extern texture<float, 1, cudaReadModeElementType> nbfp_comb_texref;
-
-/*! Texture reference for Ewald coulomb force table; bound to cu_nbparam_t.coulomb_tab */
-extern texture<float, 1, cudaReadModeElementType> coulomb_tab_texref;
-#endif /* GMX_CUDA_NB_SINGLE_COMPILATION_UNIT */
/*! Convert LJ sigma,epsilon parameters to C6,C12. */
static __forceinline__ __device__
/*! Apply potential switch, force-only version. */
static __forceinline__ __device__
void calculate_potential_switch_F(const cu_nbparam_t nbparam,
- float c6,
- float c12,
float inv_r,
float r2,
float *F_invr,
/*! Apply potential switch, force + energy version. */
static __forceinline__ __device__
void calculate_potential_switch_F_E(const cu_nbparam_t nbparam,
- float c6,
- float c12,
float inv_r,
float r2,
float *F_invr,
*E_lj *= sw;
}
+
+/*! \brief Fetch C6 grid contribution coefficients and return the product of these.
+ *
+ * Depending on what is supported, it fetches parameters either
+ * using direct load, texture objects, or texrefs.
+ */
+static __forceinline__ __device__
+float calculate_lj_ewald_c6grid(const cu_nbparam_t nbparam,
+ int typei,
+ int typej)
+{
+#if DISABLE_CUDA_TEXTURES
+ return LDG(&nbparam.nbfp_comb[2*typei]) * LDG(&nbparam.nbfp_comb[2*typej]);
+#else
+#ifdef USE_TEXOBJ
+ return tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typei) * tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typej);
+#else
+ return tex1Dfetch(nbfp_comb_texref, 2*typei) * tex1Dfetch(nbfp_comb_texref, 2*typej);
+#endif /* USE_TEXOBJ */
+#endif /* DISABLE_CUDA_TEXTURES */
+}
+
+
/*! Calculate LJ-PME grid force contribution with
* geometric combination rule.
*/
{
float c6grid, inv_r6_nm, cr2, expmcr2, poly;
-#ifdef USE_TEXOBJ
- c6grid = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typei) * tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typej);
-#else
- c6grid = tex1Dfetch(nbfp_comb_texref, 2*typei) * tex1Dfetch(nbfp_comb_texref, 2*typej);
-#endif /* USE_TEXOBJ */
+ c6grid = calculate_lj_ewald_c6grid(nbparam, typei, typej);
/* Recalculate inv_r6 without exclusion mask */
inv_r6_nm = inv_r2*inv_r2*inv_r2;
*F_invr += c6grid*(inv_r6_nm - expmcr2*(inv_r6_nm*poly + lje_coeff6_6))*inv_r2;
}
+
/*! Calculate LJ-PME grid force + energy contribution with
* geometric combination rule.
*/
{
float c6grid, inv_r6_nm, cr2, expmcr2, poly, sh_mask;
-#ifdef USE_TEXOBJ
- c6grid = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typei) * tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typej);
-#else
- c6grid = tex1Dfetch(nbfp_comb_texref, 2*typei) * tex1Dfetch(nbfp_comb_texref, 2*typej);
-#endif /* USE_TEXOBJ */
+ c6grid = calculate_lj_ewald_c6grid(nbparam, typei, typej);
/* Recalculate inv_r6 without exclusion mask */
inv_r6_nm = inv_r2*inv_r2*inv_r2;
*E_lj += c_oneSixth*c6grid*(inv_r6_nm*(1.0f - expmcr2*poly) + sh_mask);
}
+/*! Fetch per-type LJ parameters.
+ *
+ * Depending on what is supported, it fetches parameters either
+ * using direct load, texture objects, or texrefs.
+ */
+static __forceinline__ __device__
+float2 fetch_nbfp_comb_c6_c12(const cu_nbparam_t nbparam,
+ int type)
+{
+ float2 c6c12;
+#if DISABLE_CUDA_TEXTURES
+ /* Force an 8-byte fetch to save a memory instruction. */
+ float2 *nbfp_comb = (float2 *)nbparam.nbfp_comb;
+ c6c12 = LDG(&nbfp_comb[type]);
+#else
+ /* NOTE: as we always do 8-byte aligned loads, we could
+ fetch float2 here too just as above. */
+#ifdef USE_TEXOBJ
+ c6c12.x = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*type);
+ c6c12.y = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*type + 1);
+#else
+ c6c12.x = tex1Dfetch(nbfp_comb_texref, 2*type);
+ c6c12.y = tex1Dfetch(nbfp_comb_texref, 2*type + 1);
+#endif /* USE_TEXOBJ */
+#endif /* DISABLE_CUDA_TEXTURES */
+
+ return c6c12;
+}
+
+
/*! Calculate LJ-PME grid force + energy contribution (if E_lj != NULL) with
* Lorentz-Berthelot combination rule.
* We use a single F+E kernel with conditional because the performance impact
float sigma, sigma2, epsilon;
/* sigma and epsilon are scaled to give 6*C6 */
-#ifdef USE_TEXOBJ
- sigma = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typei ) + tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typej );
- epsilon = tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typei + 1) * tex1Dfetch<float>(nbparam.nbfp_comb_texobj, 2*typej + 1);
-#else
- sigma = tex1Dfetch(nbfp_comb_texref, 2*typei ) + tex1Dfetch(nbfp_comb_texref, 2*typej );
- epsilon = tex1Dfetch(nbfp_comb_texref, 2*typei + 1) * tex1Dfetch(nbfp_comb_texref, 2*typej + 1);
-#endif /* USE_TEXOBJ */
+ float2 c6c12_i = fetch_nbfp_comb_c6_c12(nbparam, typei);
+ float2 c6c12_j = fetch_nbfp_comb_c6_c12(nbparam, typej);
+
+ sigma = c6c12_i.x + c6c12_j.x;
+ epsilon = c6c12_i.y * c6c12_j.y;
+
sigma2 = sigma*sigma;
c6grid = epsilon*sigma2*sigma2*sigma2;
}
}
-/*! Interpolate Ewald coulomb force using the table through the tex_nbfp texture.
- * Original idea: from the OpenMM project
+
+/*! Fetch two consecutive values from the Ewald correction F*r table.
+ *
+ * Depending on what is supported, it fetches parameters either
+ * using direct load, texture objects, or texrefs.
*/
static __forceinline__ __device__
-float interpolate_coulomb_force_r(float r, float scale)
+float2 fetch_coulomb_force_r(const cu_nbparam_t nbparam,
+ int index)
{
- float normalized = scale * r;
- int index = (int) normalized;
- float fract2 = normalized - index;
- float fract1 = 1.0f - fract2;
+ float2 d;
+
+#if DISABLE_CUDA_TEXTURES
+ /* Can't do 8-byte fetch because some of the addresses will be misaligned. */
+ d.x = LDG(&nbparam.coulomb_tab[index]);
+ d.y = LDG(&nbparam.coulomb_tab[index + 1]);
+#else
+#ifdef USE_TEXOBJ
+ d.x = tex1Dfetch<float>(nbparam.coulomb_tab_texobj, index);
+ d.y = tex1Dfetch<float>(nbparam.coulomb_tab_texobj, index + 1);
+#else
+ d.x = tex1Dfetch(coulomb_tab_texref, index);
+ d.y = tex1Dfetch(coulomb_tab_texref, index + 1);
+#endif // USE_TEXOBJ
+#endif // DISABLE_CUDA_TEXTURES
- return fract1 * tex1Dfetch(coulomb_tab_texref, index)
- + fract2 * tex1Dfetch(coulomb_tab_texref, index + 1);
+ return d;
}
+/*! Linear interpolation using exactly two FMA operations.
+ *
+ * Implements numeric equivalent of: (1-t)*d0 + t*d1
+ * Note that CUDA does not have fnms, otherwise we'd use
+ * fma(t, d1, fnms(t, d0, d0)
+ * but input modifiers are designed for this and are fast.
+ */
+template <typename T>
+__forceinline__ __host__ __device__
+T lerp(T d0, T d1, T t)
+{
+ return fma(t, d1, fma(-t, d0, d0));
+}
+
+/*! Interpolate Ewald coulomb force correction using the F*r table.
+ */
static __forceinline__ __device__
-float interpolate_coulomb_force_r(cudaTextureObject_t texobj_coulomb_tab,
- float r, float scale)
+float interpolate_coulomb_force_r(const cu_nbparam_t nbparam,
+ float r)
{
- float normalized = scale * r;
- int index = (int) normalized;
- float fract2 = normalized - index;
- float fract1 = 1.0f - fract2;
+ float normalized = nbparam.coulomb_tab_scale * r;
+ int index = (int) normalized;
+ float fraction = normalized - index;
+
+ float2 d01 = fetch_coulomb_force_r(nbparam, index);
- return fract1 * tex1Dfetch<float>(texobj_coulomb_tab, index) +
- fract2 * tex1Dfetch<float>(texobj_coulomb_tab, index + 1);
+ return lerp(d01.x, d01.y, fraction);
}
+/*! Fetch C6 and C12 from the parameter table.
+ *
+ * Depending on what is supported, it fetches parameters either
+ * using direct load, texture objects, or texrefs.
+ */
+static __forceinline__ __device__
+void fetch_nbfp_c6_c12(float &c6,
+ float &c12,
+ const cu_nbparam_t nbparam,
+ int baseIndex)
+{
+#if DISABLE_CUDA_TEXTURES
+ /* Force an 8-byte fetch to save a memory instruction. */
+ float2 *nbfp = (float2 *)nbparam.nbfp;
+ float2 c6c12;
+ c6c12 = LDG(&nbfp[baseIndex]);
+ c6 = c6c12.x;
+ c12 = c6c12.y;
+#else
+ /* NOTE: as we always do 8-byte aligned loads, we could
+ fetch float2 here too just as above. */
+#ifdef USE_TEXOBJ
+ c6 = tex1Dfetch<float>(nbparam.nbfp_texobj, 2*baseIndex);
+ c12 = tex1Dfetch<float>(nbparam.nbfp_texobj, 2*baseIndex + 1);
+#else
+ c6 = tex1Dfetch(nbfp_texref, 2*baseIndex);
+ c12 = tex1Dfetch(nbfp_texref, 2*baseIndex + 1);
+#endif
+#endif // DISABLE_CUDA_TEXTURES
+}
+
+
/*! Calculate analytical Ewald correction term. */
static __forceinline__ __device__
float pmecorrF(float z2)
/*! Final j-force reduction; this implementation only with power of two
* array sizes and with sm >= 3.0
*/
-#if GMX_PTX_ARCH >= 300
+#if GMX_PTX_ARCH >= 300 || GMX_PTX_ARCH == 0
static __forceinline__ __device__
void reduce_force_j_warp_shfl(float3 f, float3 *fout,
- int tidxi, int aidx)
+ int tidxi, int aidx,
+ const unsigned int activemask)
{
- f.x += __shfl_down(f.x, 1);
- f.y += __shfl_up (f.y, 1);
- f.z += __shfl_down(f.z, 1);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 1);
+ f.y += gmx_shfl_up_sync (activemask, f.y, 1);
+ f.z += gmx_shfl_down_sync(activemask, f.z, 1);
if (tidxi & 1)
{
f.x = f.y;
}
- f.x += __shfl_down(f.x, 2);
- f.z += __shfl_up (f.z, 2);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 2);
+ f.z += gmx_shfl_up_sync (activemask, f.z, 2);
if (tidxi & 2)
{
f.x = f.z;
}
- f.x += __shfl_down(f.x, 4);
+ f.x += gmx_shfl_down_sync(activemask, f.x, 4);
if (tidxi < 3)
{
/*! Final i-force reduction; this implementation works only with power of two
* array sizes and with sm >= 3.0
*/
-#if GMX_PTX_ARCH >= 300
+#if GMX_PTX_ARCH >= 300 || GMX_PTX_ARCH == 0
static __forceinline__ __device__
void reduce_force_i_warp_shfl(float3 fin, float3 *fout,
float *fshift_buf, bool bCalcFshift,
- int tidxj, int aidx)
+ int tidxj, int aidx,
+ const unsigned int activemask)
{
- fin.x += __shfl_down(fin.x, c_clSize);
- fin.y += __shfl_up (fin.y, c_clSize);
- fin.z += __shfl_down(fin.z, c_clSize);
+ fin.x += gmx_shfl_down_sync(activemask, fin.x, c_clSize);
+ fin.y += gmx_shfl_up_sync (activemask, fin.y, c_clSize);
+ fin.z += gmx_shfl_down_sync(activemask, fin.z, c_clSize);
if (tidxj & 1)
{
fin.x = fin.y;
}
- fin.x += __shfl_down(fin.x, 2*c_clSize);
- fin.z += __shfl_up (fin.z, 2*c_clSize);
+ fin.x += gmx_shfl_down_sync(activemask, fin.x, 2*c_clSize);
+ fin.z += gmx_shfl_up_sync (activemask, fin.z, 2*c_clSize);
if (tidxj & 2)
{
/*! Energy reduction; this implementation works only with power of two
* array sizes and with sm >= 3.0
*/
-#if GMX_PTX_ARCH >= 300
+#if GMX_PTX_ARCH >= 300 || GMX_PTX_ARCH == 0
static __forceinline__ __device__
void reduce_energy_warp_shfl(float E_lj, float E_el,
float *e_lj, float *e_el,
- int tidx)
+ int tidx,
+ const unsigned int activemask)
{
int i, sh;
#pragma unroll 5
for (i = 0; i < 5; i++)
{
- E_lj += __shfl_down(E_lj, sh);
- E_el += __shfl_down(E_el, sh);
+ E_lj += gmx_shfl_down_sync(activemask, E_lj, sh);
+ E_el += gmx_shfl_down_sync(activemask, E_el, sh);
sh += sh;
}
"Checkpointing is merely intended for plain continuation of runs. "
"For safety reasons you must specify all file names (e.g. with -deffnm), "
"and all these files must match the names used in the run prior to checkpointing "
- "since we will append to them by default. If the files are not available, you "
+ "since we will append to them by default. If you used -deffnm and the files listed above as not "
+ "present are in fact present, try explicitly specifying them in respective mdrun options. "
+ "If the files are not available, you "
"can add the -noappend flag to mdrun and write separate new parts. "
"For mere concatenation of files, you should use the gmx trjcat tool instead.",
nfiles-nexist, nfiles);
*/
strcpy(suf_up, part_suffix);
upstring(suf_up);
- *bAddPart = (strstr(fn, part_suffix) != NULL ||
- strstr(fn, suf_up) != NULL);
+ *bAddPart = (strstr(fn, part_suffix) != nullptr ||
+ strstr(fn, suf_up) != nullptr);
}
sfree(outputfiles);
}
else
{
- fpmulti = NULL;
+ fpmulti = nullptr;
}
check_multi_int(fpmulti, cr->ms, sim_part, "simulation part", TRUE);
}
// The result will always be correct since we mask the result with m, but
// for debug builds we also want to make sure not to generate FP exceptions
#ifndef NDEBUG
- x.simdInternal_ = vbslq_f64(m.simdInternal_, x.simdInternal_, vdupq_n_f64(1.0f));
+ x.simdInternal_ = vbslq_f64(m.simdInternal_, x.simdInternal_, vdupq_n_f64(1.0));
#endif
return {
float64x2_t(vandq_u64(uint64x2_t(vrsqrteq_f64(x.simdInternal_)), m.simdInternal_))
// The result will always be correct since we mask the result with m, but
// for debug builds we also want to make sure not to generate FP exceptions
#ifndef NDEBUG
- x.simdInternal_ = vbslq_f64(m.simdInternal_, x.simdInternal_, vdupq_n_f64(1.0f));
+ x.simdInternal_ = vbslq_f64(m.simdInternal_, x.simdInternal_, vdupq_n_f64(1.0));
#endif
return {
float64x2_t(vandq_u64(uint64x2_t(vrecpeq_f64(x.simdInternal_)), m.simdInternal_))
selectByMask(SimdDouble a, SimdDBool m)
{
return {
- float64x2_t(vandq_u64(uint64x2_t(a.simdInternal_), m.simdInternal_))
+ float64x2_t(vandq_u64(uint64x2_t(a.simdInternal_), m.simdInternal_))
};
}
selectByNotMask(SimdDouble a, SimdDBool m)
{
return {
- float64x2_t(vbicq_u64(uint64x2_t(a.simdInternal_), m.simdInternal_))
+ float64x2_t(vbicq_u64(uint64x2_t(a.simdInternal_), m.simdInternal_))
};
}
blend(SimdDouble a, SimdDouble b, SimdDBool sel)
{
return {
- vbslq_f64(sel.simdInternal_, b.simdInternal_, a.simdInternal_)
+ vbslq_f64(sel.simdInternal_, b.simdInternal_, a.simdInternal_)
};
}
operator<<(SimdDInt32 a, int n)
{
return {
- vshl_n_s32(a.simdInternal_, n)
- vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? 32 : n))
++ vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? 32 : n))
};
}
operator>>(SimdDInt32 a, int n)
{
return {
- vshr_n_s32(a.simdInternal_, n)
- vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? -32 : -n))
++ vshl_s32(a.simdInternal_, vdup_n_s32(n >= 32 ? -32 : -n))
};
}
operator&(SimdDInt32 a, SimdDInt32 b)
{
return {
- vand_s32(a.simdInternal_, b.simdInternal_)
+ vand_s32(a.simdInternal_, b.simdInternal_)
};
}
andNot(SimdDInt32 a, SimdDInt32 b)
{
return {
- vbic_s32(b.simdInternal_, a.simdInternal_)
+ vbic_s32(b.simdInternal_, a.simdInternal_)
};
}
operator|(SimdDInt32 a, SimdDInt32 b)
{
return {
- vorr_s32(a.simdInternal_, b.simdInternal_)
+ vorr_s32(a.simdInternal_, b.simdInternal_)
};
}
operator^(SimdDInt32 a, SimdDInt32 b)
{
return {
- veor_s32(a.simdInternal_, b.simdInternal_)
+ veor_s32(a.simdInternal_, b.simdInternal_)
};
}
operator+(SimdDInt32 a, SimdDInt32 b)
{
return {
- vadd_s32(a.simdInternal_, b.simdInternal_)
+ vadd_s32(a.simdInternal_, b.simdInternal_)
};
}
operator-(SimdDInt32 a, SimdDInt32 b)
{
return {
- vsub_s32(a.simdInternal_, b.simdInternal_)
+ vsub_s32(a.simdInternal_, b.simdInternal_)
};
}
operator*(SimdDInt32 a, SimdDInt32 b)
{
return {
- vmul_s32(a.simdInternal_, b.simdInternal_)
+ vmul_s32(a.simdInternal_, b.simdInternal_)
};
}
operator==(SimdDInt32 a, SimdDInt32 b)
{
return {
- vceq_s32(a.simdInternal_, b.simdInternal_)
+ vceq_s32(a.simdInternal_, b.simdInternal_)
};
}
testBits(SimdDInt32 a)
{
return {
- vtst_s32( a.simdInternal_, a.simdInternal_)
+ vtst_s32( a.simdInternal_, a.simdInternal_)
};
}
operator<(SimdDInt32 a, SimdDInt32 b)
{
return {
- vclt_s32(a.simdInternal_, b.simdInternal_)
+ vclt_s32(a.simdInternal_, b.simdInternal_)
};
}
operator&&(SimdDIBool a, SimdDIBool b)
{
return {
- vand_u32(a.simdInternal_, b.simdInternal_)
+ vand_u32(a.simdInternal_, b.simdInternal_)
};
}
operator||(SimdDIBool a, SimdDIBool b)
{
return {
- vorr_u32(a.simdInternal_, b.simdInternal_)
+ vorr_u32(a.simdInternal_, b.simdInternal_)
};
}
selectByMask(SimdDInt32 a, SimdDIBool m)
{
return {
- vand_s32(a.simdInternal_, vreinterpret_s32_u32(m.simdInternal_))
+ vand_s32(a.simdInternal_, vreinterpret_s32_u32(m.simdInternal_))
};
}
selectByNotMask(SimdDInt32 a, SimdDIBool m)
{
return {
- vbic_s32(a.simdInternal_, vreinterpret_s32_u32(m.simdInternal_))
+ vbic_s32(a.simdInternal_, vreinterpret_s32_u32(m.simdInternal_))
};
}
blend(SimdDInt32 a, SimdDInt32 b, SimdDIBool sel)
{
return {
- vbsl_s32(sel.simdInternal_, b.simdInternal_, a.simdInternal_)
+ vbsl_s32(sel.simdInternal_, b.simdInternal_, a.simdInternal_)
};
}
cvtR2I(SimdDouble a)
{
return {
- vmovn_s64(vcvtnq_s64_f64(a.simdInternal_))
+ vmovn_s64(vcvtnq_s64_f64(a.simdInternal_))
};
}
cvttR2I(SimdDouble a)
{
return {
- vmovn_s64(vcvtq_s64_f64(a.simdInternal_))
+ vmovn_s64(vcvtq_s64_f64(a.simdInternal_))
};
}
cvtI2R(SimdDInt32 a)
{
return {
- vcvtq_f64_s64(vmovl_s32(a.simdInternal_))
+ vcvtq_f64_s64(vmovl_s32(a.simdInternal_))
};
}
cvtB2IB(SimdDBool a)
{
return {
- vqmovn_u64(a.simdInternal_)
+ vqmovn_u64(a.simdInternal_)
};
}
cvtIB2B(SimdDIBool a)
{
return {
- vorrq_u64(vmovl_u32(a.simdInternal_), vshlq_n_u64(vmovl_u32(a.simdInternal_), 32))
+ vorrq_u64(vmovl_u32(a.simdInternal_), vshlq_n_u64(vmovl_u32(a.simdInternal_), 32))
};
}
cvtDD2F(SimdDouble d0, SimdDouble d1)
{
return {
- vcvt_high_f32_f64(vcvt_f32_f64(d0.simdInternal_), d1.simdInternal_)
+ vcvt_high_f32_f64(vcvt_f32_f64(d0.simdInternal_), d1.simdInternal_)
};
}
#include <stdlib.h>
#include <algorithm>
+#include <memory>
#include "thread_mpi/threads.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme-load-balancing.h"
#include "gromacs/fileio/trxio.h"
-#include "gromacs/gmxlib/md_logging.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/gpu_utils/gpu_utils.h"
#include "gromacs/mdlib/mdebin.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/mdsetup.h"
#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
#include "gromacs/mdlib/ns.h"
#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
}
}
-static void reset_all_counters(FILE *fplog, t_commrec *cr,
+static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
gmx_int64_t step,
gmx_int64_t *step_rel, t_inputrec *ir,
gmx_wallcycle_t wcycle, t_nrnb *nrnb,
char sbuf[STEPSTRSIZE];
/* Reset all the counters related to performance over the run */
- md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
- gmx_step_str(step, sbuf));
+ GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+ "step %s: resetting all time and cycle counters",
+ gmx_step_str(step, sbuf));
if (use_GPU(nbv))
{
}
/*! \libinternal
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+ \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t *oenv, gmx_bool bVerbose,
int nstglobalcomm,
gmx_vsite_t *vsite, gmx_constr_t constr,
int stepout,
+ gmx::IMDOutputProvider *outputProvider,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
real cpt_period, real max_hours,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
*/
-double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
+double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+ int nfile, const t_filenm fnm[],
const gmx_output_env_t *oenv, gmx_bool bVerbose,
int nstglobalcomm,
gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout, t_inputrec *ir,
+ int stepout, gmx::IMDOutputProvider *outputProvider,
+ t_inputrec *ir,
gmx_mtop_t *top_global,
t_fcdata *fcd,
t_state *state_global,
+ ObservablesHistory *observablesHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed, t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ t_forcerec *fr,
+ const ReplicaExchangeParameters &replExParams,
gmx_membed_t *membed,
real cpt_period, real max_hours,
int imdport,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
{
- gmx_mdoutf_t outf = NULL;
+ gmx_mdoutf_t outf = nullptr;
gmx_int64_t step, step_rel;
double elapsed_time;
double t, t0, lam0[efptNR];
t_trxstatus *status;
rvec mu_tot;
t_vcm *vcm;
- matrix pcoupl_mu, M;
+ matrix parrinellorahmanMu, M;
t_trxframe rerun_fr;
- gmx_repl_ex_t repl_ex = NULL;
+ gmx_repl_ex_t repl_ex = nullptr;
int nchkpt = 1;
gmx_localtop_t *top;
- t_mdebin *mdebin = NULL;
- t_state *state = NULL;
+ t_mdebin *mdebin = nullptr;
gmx_enerdata_t *enerd;
- rvec *f = NULL;
+ PaddedRVecVector f {};
gmx_global_stat_t gstat;
- gmx_update_t *upd = NULL;
- t_graph *graph = NULL;
+ gmx_update_t *upd = nullptr;
+ t_graph *graph = nullptr;
gmx_groups_t *groups;
gmx_ekindata_t *ekind;
gmx_shellfc_t *shellfc;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bTemp, bPres, bTrotter;
real dvdl_constr;
- rvec *cbuf = NULL;
+ rvec *cbuf = nullptr;
int cbuf_nalloc = 0;
matrix lastbox;
int lamnew = 0;
/* PME load balancing data for GPU kernels */
- pme_load_balancing_t *pme_loadbal = NULL;
+ pme_load_balancing_t *pme_loadbal = nullptr;
gmx_bool bPMETune = FALSE;
gmx_bool bPMETunePrinting = FALSE;
nstglobalcomm = 1;
}
- nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
+ nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
bGStatEveryStep = (nstglobalcomm == 1);
if (bRerunMD)
}
groups = &top_global->groups;
+ gmx_edsam *ed = nullptr;
+ if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
+ {
+ /* Initialize essential dynamics sampling */
+ ed = init_edsam(opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
+ top_global, ir, cr, constr,
+ as_rvec_array(state_global->x.data()),
+ state_global->box, observablesHistory,
+ oenv, Flags & MD_APPENDFILES);
+ }
+
if (ir->eSwapCoords != eswapNO)
{
/* Initialize ion swapping code */
init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
- top_global, state_global->x, state_global->box, &state_global->swapstate, cr, oenv, Flags);
+ top_global, as_rvec_array(state_global->x.data()), state_global->box, observablesHistory, cr, oenv, Flags);
}
/* Initial values */
- init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
+ init_md(fplog, cr, outputProvider, ir, oenv, &t, &t0, state_global->lambda,
&(state_global->fep_state), lam0,
nrnb, top_global, &upd,
nfile, fnm, &outf, &mdebin,
snew(enerd, 1);
init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
enerd);
- if (DOMAINDECOMP(cr))
- {
- f = NULL;
- }
- else
- {
- snew(f, top_global->natoms);
- }
/* Kinetic energy data */
snew(ekind, 1);
}
}
+ std::unique_ptr<t_state> stateInstance;
+ t_state * state;
+
if (DOMAINDECOMP(cr))
{
top = dd_init_local_top(top_global);
- snew(state, 1);
+ stateInstance = std::unique_ptr<t_state>(new t_state);
+ state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
}
else
{
- top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
-
- state = serial_init_local_state(state_global);
-
- atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
-
- if (vsite)
- {
- set_vsite_top(vsite, top, mdatoms, cr);
- }
-
- if (ir->ePBC != epbcNONE && !fr->bMolPBC)
- {
- graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
- }
-
- if (shellfc)
- {
- make_local_shells(cr, mdatoms, shellfc);
- }
+ state_change_natoms(state_global, state_global->natoms);
+ /* We need to allocate one element extra, since we might use
+ * (unaligned) 4-wide SIMD loads to access rvec entries.
+ */
+ f.resize(state_global->natoms + 1);
+ /* Copy the pointer to the global state */
+ state = state_global;
- setup_bonded_threading(fr, &top->idef);
+ snew(top, 1);
+ mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
+ &graph, mdatoms, vsite, shellfc);
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
/* Set up interactive MD (IMD) */
- init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+ init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
nfile, fnm, oenv, imdport, Flags);
if (DOMAINDECOMP(cr))
state_global, top_global, ir,
state, &f, mdatoms, top, fr,
vsite, constr,
- nrnb, NULL, FALSE);
+ nrnb, nullptr, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
update_mdatoms(mdatoms, state->lambda[efptMASS]);
if (ir->bExpanded)
{
- init_expanded_ensemble(startingFromCheckpoint, ir, &state->dfhist);
+ init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
}
if (MASTER(cr))
/* Update mdebin with energy history if appending to output files */
if (Flags & MD_APPENDFILES)
{
- restore_energyhistory_from_state(mdebin, state_global->enerhist);
+ restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
}
- else
+ else if (observablesHistory->energyHistory.get() != nullptr)
{
- /* We might have read an energy history from checkpoint,
- * free the allocated memory and reset the counts.
+ /* We might have read an energy history from checkpoint.
+ * As we are not appending, we want to restart the statistics.
+ * Free the allocated memory and reset the counts.
*/
- done_energyhistory(state_global->enerhist);
- init_energyhistory(state_global->enerhist);
+ observablesHistory->energyHistory = {};
}
}
+ if (observablesHistory->energyHistory.get() == nullptr)
+ {
+ observablesHistory->energyHistory = std::unique_ptr<energyhistory_t>(new energyhistory_t {});
+ }
/* Set the initial energy history in state by updating once */
- update_energyhistory(state_global->enerhist, mdebin);
+ update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
}
/* Initialize constraints */
set_constraints(constr, top, ir, mdatoms, cr);
}
- if (repl_ex_nst > 0 && MASTER(cr))
+ const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
+ if (useReplicaExchange && MASTER(cr))
{
repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
- repl_ex_nst, repl_ex_nex, repl_ex_seed);
+ replExParams);
}
- /* PME tuning is only supported with PME for Coulomb. Is is not supported
- * with only LJ PME, or for reruns.
- */
+ /* PME tuning is only supported in the Verlet scheme, with PME for
+ * Coulomb. It is not supported with only LJ PME, or for
+ * reruns. */
bPMETune = ((Flags & MD_TUNEPME) && EEL_PME(fr->eeltype) && !bRerunMD &&
- !(Flags & MD_REPRODUCIBLE));
+ !(Flags & MD_REPRODUCIBLE) && ir->cutoff_scheme != ecutsGROUP);
if (bPMETune)
{
- pme_loadbal_init(&pme_loadbal, cr, fplog, ir, state->box,
- fr->ic, fr->pmedata, use_GPU(fr->nbv),
+ pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
+ fr->ic, fr->nbv->listParams.get(), fr->pmedata, use_GPU(fr->nbv),
&bPMETunePrinting);
}
if (!ir->bContinuation && !bRerunMD)
{
- if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+ if (state->flags & (1 << estV))
{
- /* Set the velocities of frozen particles to zero */
+ /* Set the velocities of vsites, shells and frozen atoms to zero */
for (i = 0; i < mdatoms->homenr; i++)
{
- for (m = 0; m < DIM; m++)
+ if (mdatoms->ptype[i] == eptVSite ||
+ mdatoms->ptype[i] == eptShell)
+ {
+ clear_rvec(state->v[i]);
+ }
+ else if (mdatoms->cFREEZE)
{
- if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ for (m = 0; m < DIM; m++)
{
- state->v[i][m] = 0;
+ if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ {
+ state->v[i][m] = 0;
+ }
}
}
}
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, state->x, ir->delta_t, NULL,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
{
nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
}
- if (repl_ex_nst > 0)
+ if (useReplicaExchange)
{
- nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+ nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
}
}
bSumEkinhOld = FALSE;
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
| (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
perhaps loses some logic?*/
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
}
// checkpoints and stop conditions to act on the same step, so
// the propagation of such signals must take place between
// simulations, not just within simulations.
- bool checkpointIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
- bool stopConditionIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+ bool checkpointIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
+ bool stopConditionIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
bool resetCountersIsLocal = true;
signals[eglsCHKPT] = SimulationSignal(checkpointIsLocal);
signals[eglsSTOPCOND] = SimulationSignal(stopConditionIsLocal);
signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
}
+ DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
+ DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
+
step = ir->init_step;
step_rel = 0;
// TODO extract this to new multi-simulation module
- if (MASTER(cr) && MULTISIM(cr) && (repl_ex_nst <= 0 ))
+ if (MASTER(cr) && MULTISIM(cr) && !useReplicaExchange)
{
if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
{
- md_print_info(cr, fplog,
- "Note: The number of steps is not consistent across multi simulations,\n"
- "but we are proceeding anyway!\n");
+ GMX_LOG(mdlog.warning).appendText(
+ "Note: The number of steps is not consistent across multi simulations,\n"
+ "but we are proceeding anyway!");
}
if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
{
- md_print_info(cr, fplog,
- "Note: The initial step is not consistent across multi simulations,\n"
- "but we are proceeding anyway!\n");
+ GMX_LOG(mdlog.warning).appendText(
+ "Note: The initial step is not consistent across multi simulations,\n"
+ "but we are proceeding anyway!");
}
}
{
/* PME grid + cut-off optimization with GPUs or PME nodes */
pme_loadbal_do(pme_loadbal, cr,
- (bVerbose && MASTER(cr)) ? stderr : NULL,
- fplog,
+ (bVerbose && MASTER(cr)) ? stderr : nullptr,
+ fplog, mdlog,
ir, fr, state,
wcycle,
step, step_rel,
&& (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
}
- bDoReplEx = ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
- do_per_step(step, repl_ex_nst));
+ bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
+ do_per_step(step, replExParams.exchangeInterval));
if (bSimAnn)
{
/* Following is necessary because the graph may get out of sync
* with the coordinates if we only have every N'th coordinate set
*/
- mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
- shift_self(graph, state->box, state->x);
+ mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
}
}
nrnb, wcycle,
do_verbose && !bPMETunePrinting);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
}
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
relax_shell_flexcon(fplog, cr, bVerbose, step,
ir, bNS, force_flags, top,
constr, enerd, fcd,
- state, f, force_vir, mdatoms,
+ state, &f, force_vir, mdatoms,
nrnb, wcycle, graph, groups,
shellfc, fr, bBornRadii, t, mu_tot,
- vsite, mdoutf_get_fp_field(outf));
+ vsite,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
}
else
{
* Check comments in sim_util.c
*/
do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
+ state->box, &state->x, &state->hist,
+ &f, force_vir, mdatoms, enerd, fcd,
state->lambda, graph,
- fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
- (bNS ? GMX_FORCE_NS : 0) | force_flags);
+ fr, vsite, mu_tot, t, ed, bBornRadii,
+ (bNS ? GMX_FORCE_NS : 0) | force_flags,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
}
if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
/* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
{
- rvec *vbuf = NULL;
+ rvec *vbuf = nullptr;
wallcycle_start(wcycle, ewcUPDATE);
if (ir->eI == eiVV && bInitStep)
* so that the input is actually the initial step.
*/
snew(vbuf, state->natoms);
- copy_rvecn(state->v, vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
}
else
{
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY1,
cr, constr);
if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
wallcycle_stop(wcycle, ewcUPDATE);
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
{
/* Need to unshift here if a do_force has been
called in the previous step */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
{
/* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
- enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
enerd->term[F_EKIN] = trace(ekind->ekin);
}
}
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
}
/* if it's the initial step, we performed this first step just to get the constraint virial */
if (ir->eI == eiVV && bInitStep)
{
- copy_rvecn(vbuf, state->v, 0, state->natoms);
+ copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
sfree(vbuf);
}
wallcycle_stop(wcycle, ewcUPDATE);
/* compute the conserved quantity */
if (EI_VV(ir->eI))
{
- saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
+ saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
if (ir->eI == eiVV)
{
last_ekin = enerd->term[F_EKIN];
statistics, but if performing simulated tempering, we
do update the velocities and the tau_t. */
- lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms);
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
- copy_df_history(&state_global->dfhist, &state->dfhist);
+ copy_df_history(state_global->dfhist, state->dfhist);
}
/* Now we have the energies and forces corresponding to the
* the update.
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
- ir, state, state_global, top_global, fr,
- outf, mdebin, ekind, f,
+ ir, state, state_global, observablesHistory,
+ top_global, fr,
+ outf, mdebin, ekind, &f,
&nchkpt,
bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
- bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
else
{
update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
- update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+ update_pcouple_before_coordinates(fplog, step, ir, state,
+ parrinellorahmanMu, M,
+ bInitStep);
}
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY2,
cr, constr);
}
cbuf_nalloc = state->natoms;
srenew(cbuf, cbuf_nalloc);
}
- copy_rvecn(state->x, cbuf, 0, state->natoms);
+ copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC, graph, f,
+ fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, lastbox,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
(bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again again */
- copy_rvecn(cbuf, state->x, 0, state->natoms);
+ copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
/* are the small terms in the shake_vir here due
* to numerical errors, or are they important
* physically? I'm thinking they are just errors, but not completely sure.
* For now, will call without actually constraining, constr=NULL*/
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
- cr, nrnb, wcycle, upd, NULL,
+ cr, nrnb, wcycle, upd, nullptr,
FALSE, bCalcVir);
}
if (EI_VV(ir->eI))
else if (graph)
{
/* Need to unshift here */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- if (vsite != NULL)
+ if (vsite != nullptr)
{
wallcycle_start(wcycle, ewcVSITECONSTR);
- if (graph != NULL)
+ if (graph != nullptr)
{
- shift_self(graph, state->box, state->x);
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
- if (graph != NULL)
+ if (graph != nullptr)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
Currently done every step so that dhdl is correct in the .edr */
sum_dhdl(enerd, state->lambda, ir->fepvals);
}
- update_box(fplog, step, ir, mdatoms, state, f,
- pcoupl_mu, nrnb, upd);
+
+ update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
+ pres, force_vir, shake_vir,
+ parrinellorahmanMu,
+ state, nrnb, upd);
/* ################# END UPDATE STEP 2 ################# */
/* #### We now have r(t+dt) and v(t+dt/2) ############# */
bSumEkinhOld = TRUE;
}
- /* ######### BEGIN PREPARING EDR OUTPUT ########### */
-
- /* use the directly determined last velocity, not actually the averaged half steps */
- if (bTrotter && ir->eI == eiVV)
+ if (bCalcEner)
{
- enerd->term[F_EKIN] = last_ekin;
- }
- enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+ /* ######### BEGIN PREPARING EDR OUTPUT ########### */
- if (EI_VV(ir->eI))
- {
- enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
- }
- else
- {
- enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
+ /* use the directly determined last velocity, not actually the averaged half steps */
+ if (bTrotter && ir->eI == eiVV)
+ {
+ enerd->term[F_EKIN] = last_ekin;
+ }
+ enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+
+ if (integratorHasConservedEnergyQuantity(ir))
+ {
+ if (EI_VV(ir->eI))
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
+ }
+ else
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
+ }
+ }
+ /* ######### END PREPARING EDR OUTPUT ########### */
}
- /* ######### END PREPARING EDR OUTPUT ########### */
/* Output stuff */
if (MASTER(cr))
if (fplog && do_log && bDoExpanded)
{
/* only needed if doing expanded ensemble */
- PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
- &state_global->dfhist, state->fep_state, ir->nstlog, step);
+ PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
+ state_global->dfhist, state->fep_state, ir->nstlog, step);
}
if (bCalcEner)
{
gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
gmx_bool do_or = do_per_step(step, ir->nstorireout);
- print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
+ print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
step, t,
eprNORMAL, mdebin, fcd, groups, &(ir->opts));
do_per_step(step, ir->swap->nstswap))
{
bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
- bRerunMD ? rerun_fr.x : state->x,
+ bRerunMD ? rerun_fr.x : as_rvec_array(state->x.data()),
bRerunMD ? rerun_fr.box : state->box,
- top_global, MASTER(cr) && bVerbose, bRerunMD);
+ MASTER(cr) && bVerbose, bRerunMD);
if (bNeedRepartition && DOMAINDECOMP(cr))
{
vsite, constr,
nrnb, wcycle, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
bFirstStep = FALSE;
/* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
- if ( (membed != NULL) && (!bLastStep) )
+ if ( (membed != nullptr) && (!bLastStep) )
{
- rescale_membed(step_rel, membed, state_global->x);
+ rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
}
if (bRerunMD)
"resetting counters later in the run, e.g. with gmx "
"mdrun -resetstep.", step);
}
- reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
- use_GPU(fr->nbv) ? fr->nbv : NULL);
+ reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
+ use_GPU(fr->nbv) ? fr->nbv : nullptr);
wcycle_set_reset_counters(wcycle, -1);
if (!(cr->duty & DUTY_PME))
{
if (bRerunMD && MASTER(cr))
{
- close_trj(status);
+ close_trx(status);
}
if (!(cr->duty & DUTY_PME))
if (bPMETune)
{
- pme_loadbal_done(pme_loadbal, cr, fplog, use_GPU(fr->nbv));
+ pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
}
done_shellfc(fplog, shellfc, step_rel);
- if (repl_ex_nst > 0 && MASTER(cr))
+ if (useReplicaExchange && MASTER(cr))
{
print_replica_exchange_statistics(fplog, repl_ex);
}
finish_swapcoords(ir->swap);
}
+ /* Do essential dynamics cleanup if needed. Close .edo file */
+ done_ed(&ed);
+
/* IMD cleanup, if bIMD is TRUE. */
IMD_finalize(ir->bIMD, ir->imd);
#include <stdio.h>
#include <string.h>
+#include <cstdlib>
+
#include "gromacs/commandline/filenm.h"
#include "gromacs/commandline/pargs.h"
-#include "gromacs/fileio/readinp.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/mdlib/main.h"
#include "gromacs/mdlib/mdrun.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/utility/arraysize.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
#include "mdrun_main.h"
+#include "repl_ex.h"
#include "runner.h"
/*! \brief Return whether either of the command-line parameters that
//! Implements C-style main function for mdrun
int gmx_mdrun(int argc, char *argv[])
+{
+ gmx::Mdrunner runner;
+ return runner.mainFunction(argc, argv);
+}
+
+namespace gmx
+{
+
+int Mdrunner::mainFunction(int argc, char *argv[])
{
const char *desc[] = {
"[THISMODULE] is the main computational chemistry engine",
"The option [TT]-pforce[tt] is useful when you suspect a simulation",
"crashes due to too large forces. With this option coordinates and",
"forces of atoms with a force larger than a certain value will",
- "be printed to stderr.",
+ "be printed to stderr. It will also terminate the run when non-finite",
+ "forces are present.",
"[PAR]",
"Checkpoints containing the complete state of the system are written",
"at regular intervals (option [TT]-cpt[tt]) to the file [TT]-cpo[tt],",
"[PAR]",
"When [TT]mdrun[tt] is started with MPI, it does not run niced by default."
};
- t_commrec *cr;
- t_filenm fnm[] = {
- { efTPR, NULL, NULL, ffREAD },
- { efTRN, "-o", NULL, ffWRITE },
- { efCOMPRESSED, "-x", NULL, ffOPTWR },
- { efCPT, "-cpi", NULL, ffOPTRD | ffALLOW_MISSING },
- { efCPT, "-cpo", NULL, ffOPTWR },
- { efSTO, "-c", "confout", ffWRITE },
- { efEDR, "-e", "ener", ffWRITE },
- { efLOG, "-g", "md", ffWRITE },
- { efXVG, "-dhdl", "dhdl", ffOPTWR },
- { efXVG, "-field", "field", ffOPTWR },
- { efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tablep", "tablep", ffOPTRD },
- { efXVG, "-tableb", "table", ffOPTRDMULT },
- { efTRX, "-rerun", "rerun", ffOPTRD },
- { efXVG, "-tpi", "tpi", ffOPTWR },
- { efXVG, "-tpid", "tpidist", ffOPTWR },
- { efEDI, "-ei", "sam", ffOPTRD },
- { efXVG, "-eo", "edsam", ffOPTWR },
- { efXVG, "-devout", "deviatie", ffOPTWR },
- { efXVG, "-runav", "runaver", ffOPTWR },
- { efXVG, "-px", "pullx", ffOPTWR },
- { efXVG, "-pf", "pullf", ffOPTWR },
- { efXVG, "-ro", "rotation", ffOPTWR },
- { efLOG, "-ra", "rotangles", ffOPTWR },
- { efLOG, "-rs", "rotslabs", ffOPTWR },
- { efLOG, "-rt", "rottorque", ffOPTWR },
- { efMTX, "-mtx", "nm", ffOPTWR },
- { efRND, "-multidir", NULL, ffOPTRDMULT},
- { efDAT, "-membed", "membed", ffOPTRD },
- { efTOP, "-mp", "membed", ffOPTRD },
- { efNDX, "-mn", "membed", ffOPTRD },
- { efXVG, "-if", "imdforces", ffOPTWR },
- { efXVG, "-swap", "swapions", ffOPTWR }
- };
- const int NFILE = asize(fnm);
- /* Command line options ! */
+ /* Command line option parameters, with their default values */
gmx_bool bDDBondCheck = TRUE;
gmx_bool bDDBondComm = TRUE;
gmx_bool bTunePME = TRUE;
- gmx_bool bVerbose = FALSE;
gmx_bool bRerunVSite = FALSE;
gmx_bool bConfout = TRUE;
gmx_bool bReproducible = FALSE;
gmx_bool bIMDterm = FALSE;
gmx_bool bIMDpull = FALSE;
- int npme = -1;
- int nstlist = 0;
- int nmultisim = 0;
- int nstglobalcomm = -1;
- int repl_ex_nst = 0;
- int repl_ex_seed = -1;
- int repl_ex_nex = 0;
- int nstepout = 100;
- int resetstep = -1;
- gmx_int64_t nsteps = -2; /* the value -2 means that the mdp option will be used */
- int imdport = 8888; /* can be almost anything, 8888 is easy to remember */
-
- rvec realddxyz = {0, 0, 0};
- const char *ddrank_opt[ddrankorderNR+1] =
- { NULL, "interleave", "pp_pme", "cartesian", NULL };
- const char *dddlb_opt[] =
- { NULL, "auto", "no", "yes", NULL };
- const char *thread_aff_opt[threadaffNR+1] =
- { NULL, "auto", "on", "off", NULL };
- const char *nbpu_opt[] =
- { NULL, "auto", "cpu", "gpu", "gpu_cpu", NULL };
- real rdd = 0.0, rconstr = 0.0, dlb_scale = 0.8, pforce = -1;
- char *ddcsx = NULL, *ddcsy = NULL, *ddcsz = NULL;
- real cpt_period = 15.0, max_hours = -1;
+ /* Command line options */
+ rvec realddxyz = {0, 0, 0};
+ const char *ddrank_opt_choices[ddrankorderNR+1] =
+ { nullptr, "interleave", "pp_pme", "cartesian", nullptr };
+ const char *dddlb_opt_choices[] =
+ { nullptr, "auto", "no", "yes", nullptr };
+ const char *thread_aff_opt_choices[threadaffNR+1] =
+ { nullptr, "auto", "on", "off", nullptr };
+ const char *nbpu_opt_choices[] =
+ { nullptr, "auto", "cpu", "gpu", "gpu_cpu", nullptr };
gmx_bool bTryToAppendFiles = TRUE;
gmx_bool bKeepAndNumCPT = FALSE;
gmx_bool bResetCountersHalfWay = FALSE;
- gmx_output_env_t *oenv = NULL;
+ const char *gpuIdTaskAssignment = "";
- /* Non transparent initialization of a complex gmx_hw_opt_t struct.
- * But unfortunately we are not allowed to call a function here,
- * since declarations follow below.
- */
- gmx_hw_opt_t hw_opt = {
- 0, 0, 0, 0, threadaffSEL, 0, 0,
- { NULL, FALSE, 0, NULL }
- };
-
- t_pargs pa[] = {
+ t_pargs pa[] = {
{ "-dd", FALSE, etRVEC, {&realddxyz},
"Domain decomposition grid, 0 is optimize" },
- { "-ddorder", FALSE, etENUM, {ddrank_opt},
+ { "-ddorder", FALSE, etENUM, {ddrank_opt_choices},
"DD rank order" },
{ "-npme", FALSE, etINT, {&npme},
"Number of separate ranks to be used for PME, -1 is guess" },
"Number of OpenMP threads per MPI rank to start (0 is guess)" },
{ "-ntomp_pme", FALSE, etINT, {&hw_opt.nthreads_omp_pme},
"Number of OpenMP threads per MPI rank to start (0 is -ntomp)" },
- { "-pin", FALSE, etENUM, {thread_aff_opt},
+ { "-pin", FALSE, etENUM, {thread_aff_opt_choices},
"Whether mdrun should try to set thread affinities" },
{ "-pinoffset", FALSE, etINT, {&hw_opt.core_pinning_offset},
"The lowest logical core number to which mdrun should pin the first thread" },
{ "-pinstride", FALSE, etINT, {&hw_opt.core_pinning_stride},
"Pinning distance in logical cores for threads, use 0 to minimize the number of threads per physical core" },
- { "-gpu_id", FALSE, etSTR, {&hw_opt.gpu_opt.gpu_id},
+ { "-gpu_id", FALSE, etSTR, {&gpuIdTaskAssignment},
"List of GPU device id-s to use, specifies the per-node PP rank to GPU mapping" },
{ "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
"Check for all bonded interactions with DD" },
"The maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
{ "-rcon", FALSE, etREAL, {&rconstr},
"Maximum distance for P-LINCS (nm), 0 is estimate" },
- { "-dlb", FALSE, etENUM, {dddlb_opt},
+ { "-dlb", FALSE, etENUM, {dddlb_opt_choices},
"Dynamic load balancing (with DD)" },
{ "-dds", FALSE, etREAL, {&dlb_scale},
"Fraction in (0,1) by whose reciprocal the initial DD cell size will be increased in order to "
"load balancing." },
{ "-gcom", FALSE, etINT, {&nstglobalcomm},
"Global communication frequency" },
- { "-nb", FALSE, etENUM, {&nbpu_opt},
+ { "-nb", FALSE, etENUM, {&nbpu_opt_choices},
"Calculate non-bonded interactions on" },
- { "-nstlist", FALSE, etINT, {&nstlist},
+ { "-nstlist", FALSE, etINT, {&nstlist_cmdline},
"Set nstlist when using a Verlet buffer tolerance (0 is guess)" },
{ "-tunepme", FALSE, etBOOL, {&bTunePME},
- "Optimize PME load between PP/PME ranks or GPU/CPU" },
+ "Optimize PME load between PP/PME ranks or GPU/CPU (only with the Verlet cut-off scheme)" },
{ "-v", FALSE, etBOOL, {&bVerbose},
"Be loud and noisy" },
{ "-pforce", FALSE, etREAL, {&pforce},
"Keep and number checkpoint files" },
{ "-append", FALSE, etBOOL, {&bTryToAppendFiles},
"Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names" },
- { "-nsteps", FALSE, etINT64, {&nsteps},
+ { "-nsteps", FALSE, etINT64, {&nsteps_cmdline},
"Run this number of steps, overrides .mdp file option (-1 means infinite, -2 means use mdp option, smaller is invalid)" },
{ "-maxh", FALSE, etREAL, {&max_hours},
"Terminate after 0.99 times this time (hours)" },
{ "-multi", FALSE, etINT, {&nmultisim},
"Do multiple simulations in parallel" },
- { "-replex", FALSE, etINT, {&repl_ex_nst},
+ { "-replex", FALSE, etINT, {&replExParams.exchangeInterval},
"Attempt replica exchange periodically with this period (steps)" },
- { "-nex", FALSE, etINT, {&repl_ex_nex},
+ { "-nex", FALSE, etINT, {&replExParams.numExchanges},
"Number of random exchanges to carry out each exchange interval (N^3 is one suggestion). -nex zero or not specified gives neighbor replica exchange." },
- { "-reseed", FALSE, etINT, {&repl_ex_seed},
+ { "-reseed", FALSE, etINT, {&replExParams.randomSeed},
"Seed for replica exchange, -1 is generate a seed" },
{ "-imdport", FALSE, etINT, {&imdport},
"HIDDENIMD listening port" },
{ "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
"HIDDENReset the cycle counters after half the number of steps or halfway [TT]-maxh[tt]" }
};
- unsigned long Flags;
- ivec ddxyz;
- int dd_rank_order;
- gmx_bool bDoAppendFiles, bStartFromCpt;
- FILE *fplog;
- int rc;
- char **multidir = NULL;
+ gmx_bool bDoAppendFiles, bStartFromCpt;
+ int rc;
+ char **multidir = nullptr;
cr = init_commrec();
}
*/
- if (!parse_common_args(&argc, argv, PCA_Flags, NFILE, fnm, asize(pa), pa,
- asize(desc), desc, 0, NULL, &oenv))
+ if (!parse_common_args(&argc, argv, PCA_Flags, nfile, fnm, asize(pa), pa,
+ asize(desc), desc, 0, nullptr, &oenv))
{
+ sfree(cr);
return 0;
}
+ // Handle the option that permits the user to select a GPU task
+ // assignment, which could be in an environment variable (so that
+ // there is a way to customize it, when using MPI in heterogeneous
+ // contexts).
+ {
+ // TODO Argument parsing can't handle std::string. We should
+ // fix that by changing the parsing, once more of the roles of
+ // handling, validating and implementing defaults for user
+ // command-line options have been seperated.
+ hw_opt.gpuIdTaskAssignment = gpuIdTaskAssignment;
+ const char *env = getenv("GMX_GPU_ID");
+ if (env != nullptr)
+ {
+ if (!hw_opt.gpuIdTaskAssignment.empty())
+ {
+ gmx_fatal(FARGS, "GMX_GPU_ID and -gpu_id can not be used at the same time");
+ }
+ hw_opt.gpuIdTaskAssignment = env;
+ }
+ }
- dd_rank_order = nenum(ddrank_opt);
+ dd_rank_order = nenum(ddrank_opt_choices);
- hw_opt.thread_affinity = nenum(thread_aff_opt);
+ hw_opt.thread_affinity = nenum(thread_aff_opt_choices);
/* now check the -multi and -multidir option */
- if (opt2bSet("-multidir", NFILE, fnm))
+ if (opt2bSet("-multidir", nfile, fnm))
{
if (nmultisim > 0)
{
gmx_fatal(FARGS, "mdrun -multi and -multidir options are mutually exclusive.");
}
- nmultisim = opt2fns(&multidir, "-multidir", NFILE, fnm);
+ nmultisim = opt2fns(&multidir, "-multidir", nfile, fnm);
}
- if (repl_ex_nst != 0 && nmultisim < 2)
+ if (replExParams.exchangeInterval != 0 && nmultisim < 2)
{
gmx_fatal(FARGS, "Need at least two replicas for replica exchange (option -multi)");
}
- if (repl_ex_nex < 0)
+ if (replExParams.numExchanges < 0)
{
gmx_fatal(FARGS, "Replica exchange number of exchanges needs to be positive");
}
if (nmultisim >= 1)
{
#if !GMX_THREAD_MPI
- init_multisystem(cr, nmultisim, multidir, NFILE, fnm);
+ init_multisystem(cr, nmultisim, multidir, nfile, fnm);
#else
gmx_fatal(FARGS, "mdrun -multi or -multidir are not supported with the thread-MPI library. "
"Please compile GROMACS with a proper external MPI library.");
#endif
}
- if (!opt2bSet("-cpi", NFILE, fnm))
+ if (!opt2bSet("-cpi", nfile, fnm))
{
// If we are not starting from a checkpoint we never allow files to be appended
// to, since that has caused a ton of strange behaviour and bugs in the past.
}
}
- handleRestart(cr, bTryToAppendFiles, NFILE, fnm, &bDoAppendFiles, &bStartFromCpt);
+ handleRestart(cr, bTryToAppendFiles, nfile, fnm, &bDoAppendFiles, &bStartFromCpt);
- Flags = opt2bSet("-rerun", NFILE, fnm) ? MD_RERUN : 0;
+ Flags = opt2bSet("-rerun", nfile, fnm) ? MD_RERUN : 0;
Flags = Flags | (bDDBondCheck ? MD_DDBONDCHECK : 0);
Flags = Flags | (bDDBondComm ? MD_DDBONDCOMM : 0);
Flags = Flags | (bTunePME ? MD_TUNEPME : 0);
there instead. */
if (MASTER(cr) && !bDoAppendFiles)
{
- gmx_log_open(ftp2fn(efLOG, NFILE, fnm), cr,
+ gmx_log_open(ftp2fn(efLOG, nfile, fnm), cr,
Flags & MD_APPENDFILES, &fplog);
}
else
{
- fplog = NULL;
+ fplog = nullptr;
}
ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
- rc = gmx::mdrunner(&hw_opt, fplog, cr, NFILE, fnm, oenv, bVerbose,
- nstglobalcomm, ddxyz, dd_rank_order, npme, rdd, rconstr,
- dddlb_opt[0], dlb_scale, ddcsx, ddcsy, ddcsz,
- nbpu_opt[0], nstlist,
- nsteps, nstepout, resetstep,
- nmultisim, repl_ex_nst, repl_ex_nex, repl_ex_seed,
- pforce, cpt_period, max_hours, imdport, Flags);
+ dddlb_opt = dddlb_opt_choices[0];
+ nbpu_opt = nbpu_opt_choices[0];
+ rc = mdrunner();
/* Log file has to be closed in mdrunner if we are appending to it
(fplog not set here) */
return rc;
}
+
+} // namespace
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff