-cmake_minimum_required(VERSION 2.6)
+cmake_minimum_required(VERSION 2.8)
# Keep CMake suitably quiet on Cygwin
set(CMAKE_LEGACY_CYGWIN_WIN32 0) # Remove when CMake >= 2.8.4 is required
+# override bugs on OS X where Cmake picks gcc (GNU) for C instead of system default cc (Clang).
+if(APPLE)
+ set(CMAKE_C_COMPILER_INIT "cc")
+endif(APPLE)
+
project(Gromacs)
include(Dart)
mark_as_advanced(DART_ROOT)
# machine with no git.
#
# NOTE: when releasing the "-dev" suffix needs to be stripped off!
-set(PROJECT_VERSION "4.5.5-dev")
+set(PROJECT_VERSION "4.6-dev")
set(CUSTOM_VERSION_STRING ""
CACHE STRING "Custom version string (if empty, use hard-coded default)")
mark_as_advanced(CUSTOM_VERSION_STRING)
set(SOVERSION 6)
# It is a bit irritating, but this has to be set separately for now!
SET(CPACK_PACKAGE_VERSION_MAJOR "4")
-SET(CPACK_PACKAGE_VERSION_MINOR "5")
-SET(CPACK_PACKAGE_VERSION_PATCH "5")
+SET(CPACK_PACKAGE_VERSION_MINOR "6")
+#SET(CPACK_PACKAGE_VERSION_PATCH "0")
# Cmake modules/macros are in a subdirectory to keep this file cleaner
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
-if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
+if(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT AND UNIX)
set(CMAKE_INSTALL_PREFIX "/usr/local/gromacs" CACHE STRING "Installation prefix (installation will need write permissions here)" FORCE)
endif()
enable_language(C)
+set(GMX_USE_RELATIVE_INSTALL_PATH OFF CACHE STRING "Use relative paths not absolute paths for cmake install. Has only an effect on cpack.")
+mark_as_advanced(GMX_USE_RELATIVE_INSTALL_PATH)
+
set(CPACK_PACKAGE_VERSION ${PROJECT_VERSION})
set(CPACK_PACKAGE_VENDOR "gromacs.org")
set(CPACK_PACKAGE_DESCRIPTION_SUMMARY "Gromacs - a toolkit for high-performance molecular simulation")
-set(CPACK_SET_DESTDIR "ON")
+if (NOT GMX_USE_RELATIVE_INSTALL_PATH)
+ set(CPACK_SET_DESTDIR "ON")
+endif()
set(CPACK_INCLUDE_TOPLEVEL_DIRECTORY 0)
set(CPACK_RESOURCE_FILE_WELCOME "${CMAKE_SOURCE_DIR}/admin/InstallWelcome.txt")
# Its GPL/LGPL, so they do not have to agree to a license for mere usage, but some installers require this...
set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_SOURCE_DIR}/admin/InstallLicense.txt")
set(CPACK_RESOURCE_FILE_README "${CMAKE_SOURCE_DIR}/admin/InstallInfo.txt")
+set(CPACK_SOURCE_IGNORE_FILES ".isreposource;.git;.gitignore")
+set(CPACK_PROJECT_CONFIG_FILE "${CMAKE_SOURCE_DIR}/CPackInit.cmake")
+SET(CPACK_SOURCE_INSTALLED_DIRECTORIES "${CMAKE_SOURCE_DIR};/;${CMAKE_BINARY_DIR}/man;man")
#must come after all cpack settings!
include(CPack)
# Fix stupid flags on Windows
########################################################################
SET(SHARED_LIBS_DEFAULT ON)
-IF( MSVC )
+IF( WIN32 AND NOT CYGWIN)
+ option(GMX_PREFER_STATIC_LIBS "When finding libraries prefer static system libraries (MT instead of MD)!" ON)
+ mark_as_advanced(GMX_PREFER_STATIC_LIBS)
+ SET(SHARED_LIBS_DEFAULT OFF) #is currently not working on Windows
+
+ IF (GMX_PREFER_STATIC_LIBS)
+ #Only setting Debug and Release flags. Others configurations current not used.
STRING(REPLACE /MD /MT CMAKE_C_FLAGS_RELEASE ${CMAKE_C_FLAGS_RELEASE})
SET(CMAKE_C_FLAGS_RELEASE ${CMAKE_C_FLAGS_RELEASE} CACHE STRING "" FORCE)
STRING(REPLACE /MD /MT CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG})
- SET(CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG} CACHE STRING "" FORCE)
- SET(SHARED_LIBS_DEFAULT OFF)
-ENDIF()
-
-IF( WIN32 AND NOT CYGWIN AND CMAKE_C_COMPILER_ID MATCHES "Intel" )
- STRING(REPLACE /GZ /RTC1 CMAKE_C_FLAGS_RELEASE ${CMAKE_C_FLAGS_RELEASE})
- SET(CMAKE_C_FLAGS_RELEASE ${CMAKE_C_FLAGS_RELEASE} CACHE STRING "" FORCE)
+ SET(CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG} CACHE STRING "" FORCE)
+ STRING(REPLACE /MD /MT CMAKE_CXX_FLAGS_RELEASE ${CMAKE_CXX_FLAGS_RELEASE})
+ SET(CMAKE_CXX_FLAGS_RELEASE ${CMAKE_CXX_FLAGS_RELEASE} CACHE STRING "" FORCE)
+ STRING(REPLACE /MD /MT CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG})
+ SET(CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG} CACHE STRING "" FORCE)
+ ENDIF()
+
+ #Workaround for cmake bug 13174. Replace deprecated options.
+ IF( CMAKE_C_COMPILER_ID MATCHES "Intel" )
STRING(REPLACE /GZ /RTC1 CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG})
- SET(CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG} CACHE STRING "" FORCE)
- SET(SHARED_LIBS_DEFAULT OFF)
+ SET(CMAKE_C_FLAGS_DEBUG ${CMAKE_C_FLAGS_DEBUG} CACHE STRING "" FORCE)
+ ENDIF()
+ IF( CMAKE_CXX_COMPILER_ID MATCHES "Intel" )
+ STRING(REPLACE /GZ /RTC1 CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG})
+ STRING(REPLACE /GX /EHsc CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG})
+ SET(CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG} CACHE STRING "" FORCE)
+
+ STRING(REPLACE /GX /EHsc CMAKE_CXX_FLAGS_RELEASE ${CMAKE_CXX_FLAGS_RELEASE})
+ SET(CMAKE_CXX_FLAGS_RELEASE ${CMAKE_CXX_FLAGS_RELEASE} CACHE STRING "" FORCE)
+ ENDIF()
ENDIF()
-set(GMX_EXTRA_LIBRARIES)
-
-
-
-######################################################################
-# compiler tests
-# these need ot be done early (before further tests).
-#####################################################################
-
-include(CheckCCompilerFlag)
-include(CheckCXXCompilerFlag)
-
-include(gmxCFlags)
-gmx_c_flags()
-
########################################################################
# User input options #
########################################################################
option(GMX_DOUBLE "Use double precision (much slower, use only if you really need it)" OFF)
option(GMX_MPI "Build a parallel (message-passing) version of GROMACS" OFF)
-option(GMX_THREADS "Build a parallel (thread-based) version of GROMACS (cannot be combined with MPI yet)" ON)
+option(GMX_THREAD_MPI "Build a thread-MPI-based multithreaded version of GROMACS (not compatible with MPI)" ON)
option(GMX_SOFTWARE_INVSQRT "Use GROMACS software 1/sqrt" ON)
mark_as_advanced(GMX_SOFTWARE_INVSQRT)
option(GMX_POWERPC_INVSQRT "Use PowerPC hardware 1/sqrt" OFF)
option(GMX_FAHCORE "Build a library with mdrun functionality" OFF)
mark_as_advanced(GMX_FAHCORE)
option(GMX_OPENMM "Accelerated execution on GPUs through the OpenMM library (rerun cmake after changing to see relevant options)" OFF)
-set(GMX_ACCELERATION "auto"
- CACHE STRING "Accelerated kernels. Pick one of: auto, none, SSE, BlueGene, Power6, ia64, altivec, fortran")
+
+include(gmxDetectAcceleration)
+if(NOT DEFINED GMX_ACCELERATION)
+ if(CMAKE_CROSSCOMPILING)
+ set(GMX_SUGGESTED_ACCELERATION "None")
+ else(CMAKE_CROSSCOMPILING)
+ gmx_detect_acceleration(GMX_SUGGESTED_ACCELERATION)
+ endif(CMAKE_CROSSCOMPILING)
+endif(NOT DEFINED GMX_ACCELERATION)
+
+set(GMX_ACCELERATION "@GMX_SUGGESTED_ACCELERATION@"
+ CACHE STRING "Accelerated kernels. Pick one of: None, SSE2, SSE4.1, AVX_128_FMA, AVX_256, BlueGene, Power6, Fortran")
set(GMX_FFT_LIBRARY "fftw3"
- CACHE STRING "FFT library choices: fftw3,fftw2,mkl,fftpack[built-in]")
+ CACHE STRING "FFT library choices: fftw3,mkl,fftpack[built-in]")
option(GMX_DISABLE_FFTW_MEASURE
"Do not optimize FFTW setups (not needed with SSE)" OFF)
mark_as_advanced(GMX_DISABLE_FFTW_MEASURE)
CACHE STRING "QM package choices: none,gaussian,mopac,gamess,orca")
option(GMX_BROKEN_CALLOC "Work around broken calloc()" OFF)
mark_as_advanced(GMX_BROKEN_CALLOC)
-option(BUILD_SHARED_LIBS "Enable shared libraries (can be problematic with MPI, Windows)" ${SHARED_LIBS_DEFAULT})
option(GMX_MPI_IN_PLACE "Enable MPI_IN_PLACE for MPIs that have it defined" ON)
mark_as_advanced(GMX_MPI_IN_PLACE)
-option(GMX_DLOPEN "Compile with dlopen, needed to read VMD supported file formats" ON)
-mark_as_advanced(GMX_DLOPEN)
+option(GMX_LOAD_PLUGINS "Compile with plugin support, needed to read VMD supported file formats" ON)
+mark_as_advanced(GMX_LOAD_PLUGINS)
-
-option(GMX_IA32_ASM "Add SSE assembly files for IA32" OFF)
-mark_as_advanced(GMX_IA32_ASM)
-option(GMX_X86_64_ASM "Add SSE assembly files for X86_64" OFF)
-mark_as_advanced(GMX_X86_64_ASM)
+option(GMX_OPENMP "Enable OpenMP-based mutithreading. " ON)
option(USE_VERSION_H "Generate development version string/information" ON)
mark_as_advanced(USE_VERSION_H)
mark_as_advanced(GMX_PREFER_STATIC_LIBS)
endif()
+
+######################################################################
+# compiler tests
+# these need ot be done early (before further tests).
+#####################################################################
+
+# cmake/Check{C,CXX}CompilerFlag.cmake are lifted from CMake git next
+# branch (proposed for v2.8.9) to be able to detect invalid options
+# with the Intel Compilers.
+# Remove these files from the source tree when a CMake version that
+# includes the features in question becomes required.
+include(CheckCCompilerFlag)
+include(CheckCXXCompilerFlag)
+
+# OpenMP check must come before other CFLAGS!
+if(GMX_OPENMP)
+ find_package(OpenMP)
+ if(OPENMP_FOUND)
+ set(GROMACS_C_FLAGS "${OpenMP_C_FLAGS} ${GROMACS_C_FLAGS}")
+ set(GROMACS_CXX_FLAGS "${OpenMP_CXX_FLAGS} ${GROMACS_CXX_FLAGS}")
+ else(OPENMP_FOUND)
+ message(WARNING
+ "Compiler not supporting OpenMP. This might hurt your performance a lot, "
+ "in particular with GPUs. Try using a different compiler (Intel is good on x86) "
+ "if you can - for now we are proceeding by turning off OpenMP.")
+ set(GMX_OPENMP OFF CACHE STRING "Compiler does not support OpenMP." FORCE)
+ endif(OPENMP_FOUND)
+endif()
+
+
+include(gmxCFlags)
+gmx_c_flags()
+
########################################################################
# Set up binary and library suffixing
########################################################################
set (GMX_LIBS_SUFFIX "${GMX_LIBS_SUFFIX}_d")
endif(GMX_DOUBLE)
if (GMX_OPENMM)
- set (GMX_BINARY_SUFFIX "-gpu")
- set (GMX_LIBS_SUFFIX "_gpu")
+ set (GMX_BINARY_SUFFIX "-openmm")
+ set (GMX_LIBS_SUFFIX "_openmm")
endif(GMX_OPENMM)
mark_as_advanced(FORCE GMX_BINARY_SUFFIX GMX_LIBS_SUFFIX)
- message(STATUS "Using default binary suffix: \"${GMX_BINARY_SUFFIX}\"")
- message(STATUS "Using default library suffix: \"${GMX_LIBS_SUFFIX}\"")
+ if (NOT SUFFIX_QUIETLY)
+ message(STATUS "Using default binary suffix: \"${GMX_BINARY_SUFFIX}\"")
+ message(STATUS "Using default library suffix: \"${GMX_LIBS_SUFFIX}\"")
+ endif (NOT SUFFIX_QUIETLY)
else(GMX_DEFAULT_SUFFIX)
mark_as_advanced(CLEAR GMX_BINARY_SUFFIX GMX_LIBS_SUFFIX)
- message(STATUS "Using manually set binary suffix: \"${GMX_BINARY_SUFFIX}\"")
- message(STATUS "Using manually set library suffix: \"${GMX_LIBS_SUFFIX}\"")
+ if (NOT SUFFIX_QUIETLY)
+ message(STATUS "Using manually set binary suffix: \"${GMX_BINARY_SUFFIX}\"")
+ message(STATUS "Using manually set library suffix: \"${GMX_LIBS_SUFFIX}\"")
+ endif (NOT SUFFIX_QUIETLY)
endif(GMX_DEFAULT_SUFFIX)
+set(SUFFIX_QUIETLY TRUE CACHE INTERNAL "")
set(PKG_CFLAGS "")
if(GMX_DOUBLE)
########################################################################
#Process MPI settings
########################################################################
-include(CheckCSourceCompiles) # for gmxTestMPI_IN_PLACE
-if(GMX_MPI)
- if(GMX_THREADS)
- #message(FATAL_ERROR "Thread-based parallelization conflicts with MPI.")
- set(GMX_THREADS OFF CACHE BOOL
- "Thread-based parallelization conflicts with MPI." FORCE)
- endif(GMX_THREADS)
- find_package(MPI)
- if(MPI_FOUND)
- set(GROMACS_C_FLAGS ${GROMACS_FLAGS} ${MPI_COMPILE_FLAGS})
- set(GROMACS_LINKER_FLAGS ${GROMACS_LINKER_FLAGS} ${MPI_LINK_FLAGS})
-
- include_directories(${MPI_INCLUDE_PATH})
- list(APPEND GMX_EXTRA_LIBRARIES ${MPI_LIBRARIES})
- if(GMX_FAHCORE)
- add_definitions( -DMPI ) #for FAHCORE
- endif(GMX_FAHCORE)
- include(gmxTestMPI_IN_PLACE)
- if (GMX_MPI_IN_PLACE)
- gmx_test_mpi_in_place(MPI_IN_PLACE_EXISTS)
- endif (GMX_MPI_IN_PLACE)
-
- # test for unsuitable versions of OpenMPI
- exec_program(${MPIEXEC}
- ARGS --version
- OUTPUT_VARIABLE MPI_TYPE
- RETURN_VALUE MPI_EXEC_RETURN)
- if(MPI_EXEC_RETURN EQUAL 0)
- if(MPI_TYPE MATCHES "Open MPI|OpenRTE")
- string(REGEX MATCH "[0-9]+\\.[0-9]*\\.?[0-9]*" MPI_VERSION ${MPI_TYPE})
- if(MPI_VERSION VERSION_LESS "1.4.1")
- MESSAGE(WARNING "
- There are known problems with OpenMPI version < 1.4.1.
- Please consider updating your OpenMPI.")
- endif(MPI_VERSION VERSION_LESS "1.4.1")
- unset(MPI_VERSION)
- else(MPI_TYPE MATCHES "Open MPI|OpenRTE")
- # This is not OpenMPI, so give the old generic warning message
- MESSAGE(WARNING "
- There are known problems with some MPI implementations:
- MVAPICH2 version <= 1.4.1
- Please consider updating your MPI if applicable.")
- endif(MPI_TYPE MATCHES "Open MPI|OpenRTE")
- unset(MPI_TYPE)
- endif(MPI_EXEC_RETURN EQUAL 0)
- else(MPI_FOUND)
- message(FATAL_ERROR "MPI support requested, but no MPI compiler found.")
- endif(MPI_FOUND)
- include(gmxTestCatamount)
- gmx_test_catamount(GMX_CRAY_XT3)
- if(GMX_CRAY_XT3)
- set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_CRAY_XT3")
- set(GMX_NO_SYSTEM 1)
- set(GMX_NO_NICE 1)
- endif(GMX_CRAY_XT3)
- set(GMX_LIB_MPI 1)
- set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_LIB_MPI")
-endif(GMX_MPI)
-
+include(gmxManageMPI)
#######################################################################
# Check for options incompatible with OpenMM build #
#######################################################################
if(GMX_OPENMM)
- cmake_minimum_required(VERSION 2.6.4)
# we'll use the built-in fft to avoid unnecessary dependencies
string(TOUPPER ${GMX_FFT_LIBRARY} GMX_FFT_LIBRARY)
if(NOT ${GMX_FFT_LIBRARY} STREQUAL "FFTPACK")
if(GMX_MPI)
message(FATAL_ERROR "The OpenMM build is not compatible with MPI!")
endif(GMX_MPI)
- if(GMX_THREADS)
- message(STATUS "Threads are not compatible with OpenMM build, disabled!")
- set(GMX_THREADS OFF CACHE BOOL
- "Threads are not compatible with OpenMM build, disabled!" FORCE)
- endif(GMX_THREADS)
+ if(GMX_THREAD_MPI)
+ message(STATUS "Thread-MPI not compatible with OpenMM, disabled!")
+ set(GMX_THREAD_MPI OFF CACHE BOOL
+ "Thread-MPI not compatible with OpenMM build, disabled!" FORCE)
+ endif(GMX_THREAD_MPI)
+ if(GMX_OPENMP)
+ message(STATUS "OpenMP multithreading not compatible with OpenMM, disabled")
+ set(GMX_OPENMP OFF CACHE BOOL
+ "OpenMP multithreading not compatible with OpenMM, disabled!" FORCE)
+ endif()
if(GMX_SOFTWARE_INVSQRT)
set(GMX_SOFTWARE_INVSQRT OFF CACHE STRING
"The OpenMM build does not need GROMACS software 1/sqrt!" FORCE)
endif()
# mark as advanced the unused variables
mark_as_advanced(FORCE GMX_ACCELERATION GMX_MPI GMX_FFT_LIBRARY
- GMX_QMMM_PROGRAM GMX_THREADS GMX_DOUBLE)
+ GMX_QMMM_PROGRAM GMX_THREAD_MPI GMX_DOUBLE)
else(GMX_OPENMM)
mark_as_advanced(CLEAR GMX_ACCELERATION GMX_MPI GMX_FFT_LIBRARY
- GMX_QMMM_PROGRAM GMX_THREADS GMX_DOUBLE)
+ GMX_QMMM_PROGRAM GMX_THREAD_MPI GMX_DOUBLE)
endif(GMX_OPENMM)
+
+
########################################################################
# Basic system tests (standard libraries, headers, functions, types) #
########################################################################
check_include_files(sys/time.h HAVE_SYS_TIME_H)
check_include_files(rpc/rpc.h HAVE_RPC_RPC_H)
check_include_files("rpc/rpc.h;rpc/xdr.h" HAVE_RPC_XDR_H)
-# SSE support
-check_include_files(xmmintrin.h HAVE_XMMINTRIN_H)
-check_include_files(emmintrin.h HAVE_EMMINTRIN_H)
-check_include_files(pmmintrin.h HAVE_PMMINTRIN_H)
-check_include_files(smmintrin.h HAVE_SMMINTRIN_H)
check_include_files(io.h HAVE_IO_H)
check_function_exists(fileno HAVE_FILENO)
check_function_exists(_commit HAVE__COMMIT)
check_function_exists(lstat HAVE_LSTAT)
+check_function_exists(sigaction HAVE_SIGACTION)
include(CheckLibraryExists)
check_library_exists(m sqrt "" HAVE_LIBM)
if(UNIX AND NOT APPLE)
if(GMX_PREFER_STATIC_LIBS)
SET(CMAKE_FIND_LIBRARY_SUFFIXES .a ${CMAKE_FIND_LIBRARY_SUFFIXES})
- if(BUILD_SHARED_LIBS)
- message(WARNING "Static libraries requested, the GROMACS libraries will also be build static (BUILD_SHARED_LIBS=OFF)")
- set(BUILD_SHARED_LIBS OFF CACHE BOOL "Enable shared libraries (can be problematic with MPI, Windows)" FORCE)
+ if(SHARED_LIBS_DEFAULT)
+ if(BUILD_SHARED_LIBS) #Warn the user about the combination. But don't overwrite the request.
+ message(WARNING "Static libraries requested, and shared Gromacs libraries requested.")
+ elseif(NOT DEFINED BUILD_SHARED_LIBS) #Change default to OFF. Don't warn if it's already off.
+ message(WARNING "Static libraries requested, the GROMACS libraries will also be build static (BUILD_SHARED_LIBS=OFF)")
+ set(SHARED_LIBS_DEFAULT OFF)
+ endif()
endif()
endif()
endif()
+option(BUILD_SHARED_LIBS "Enable shared libraries (can be problematic with MPI, Windows)" ${SHARED_LIBS_DEFAULT})
option(GMX_XML "Use libxml2 to parse xml files" ON)
if (GMX_XML)
endif(X11_FOUND)
endif(GMX_X11)
-if(GMX_THREADS)
- set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_THREADS")
+if(GMX_THREAD_MPI)
+ set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_THREAD_MPI")
include(ThreadMPI)
set(THREAD_MPI_LIB thread_mpi)
set(GMX_MPI 1)
string(TOUPPER ${GMX_FFT_LIBRARY} ${GMX_FFT_LIBRARY})
- if(${GMX_FFT_LIBRARY} STREQUAL "FFTW2")
- message(FATAL_ERROR "FFTW2 can't be used with threads. Try fftw3 or mkl.")
- endif()
-endif(GMX_THREADS)
+endif(GMX_THREAD_MPI)
if(GMX_OPENMM)
set(CUDA_BUILD_EMULATION OFF)
find_package(OpenMM)
endif(GMX_OPENMM)
+
if(APPLE)
find_library(ACCELERATE_FRAMEWORK Accelerate)
list(APPEND GMX_EXTRA_LIBRARIES ${ACCELERATE_FRAMEWORK})
set(GMX_CYGWIN 1)
endif(CYGWIN)
+if(WIN32 AND NOT CYGWIN)
+ set(GMX_NATIVE_WINDOWS 1)
+endif()
+
# only bother with finding git and using version.h if the source is a git repo
if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
if(USE_VERSION_H)
# We need at least git v1.5.3 be able to parse git's date output. If not
- # fund or the version is too small, we can't generate version information.
- find_package(Git 1.5.3)
- # this should at some point become VERSION_LESS
- if(NOT Git_FOUND OR Git_VERSION STRLESS "1.5.3")
- message("No compatible git version found, won't be able to generate proper development version information.")
- set(USE_VERSION_H OFF)
+ # found or the version is too small, we can't generate version information.
+ find_package(Git)
+
+ # Find out the git version
+ if(GIT_FOUND AND NOT GIT_VERSION)
+ execute_process(COMMAND ${GIT_EXECUTABLE} "--version"
+ OUTPUT_VARIABLE _exec_out
+ OUTPUT_STRIP_TRAILING_WHITESPACE)
+ string(REGEX REPLACE "git version (.*)" "\\1" GIT_VERSION ${_exec_out})
+ set(GIT_VERSION ${GIT_VERSION} CACHE STRING "Git version")
+ mark_as_advanced(GIT_VERSION)
+ endif()
+
+ if(NOT GIT_FOUND OR GIT_VERSION VERSION_LESS "1.5.3")
+ message("No compatible git version found, won't be able to generate proper development version information.")
+ set(USE_VERSION_H OFF)
endif()
endif()
else()
set(USE_VERSION_H OFF)
endif()
-if (GMX_DLOPEN)
- list(APPEND GMX_EXTRA_LIBRARIES ${CMAKE_DL_LIBS})
- set(PKG_DL_LIBS "-l${CMAKE_DL_LIBS}")
-else(GMX_DLOPEN)
- set(PKG_DL_LIBS)
-endif (GMX_DLOPEN)
-
########################################################################
# Generate development version info for cache
########################################################################
add_definitions( -DHAVE_CONFIG_H )
include_directories(${CMAKE_BINARY_DIR}/src)
+include_directories(${CMAKE_BINARY_DIR}/include)
include_directories(${CMAKE_SOURCE_DIR}/include)
-include(gmxCheckBuildUserTime)
-gmx_check_build_user_time(BUILD_TIME BUILD_USER BUILD_MACHINE)
+include(gmxTestInlineASM)
+gmx_test_inline_asm_gcc_x86(GMX_X86_GCC_INLINE_ASM)
+
+include(gmxSetBuildInformation)
+gmx_set_build_information()
+if(BUILD_CPU_FEATURES MATCHES "rdtscp")
+ # The timestep counter headers do not include config.h
+ add_definitions(-DHAVE_RDTSCP)
+endif(BUILD_CPU_FEATURES MATCHES "rdtscp")
include(gmxTestFloatFormat)
gmx_test_float_format(GMX_FLOAT_FORMAT_IEEE754
check_function_exists(fseeko HAVE_FSEEKO)
include(gmxTestSignal)
-gmx_test_retsigtype(RETSIGTYPE)
gmx_test_sigusr1(HAVE_SIGUSR1)
include(gmxTestInline)
gmx_test__isfinite(HAVE__ISFINITE)
gmx_test__finite(HAVE__FINITE)
-include(gmxTestInlineASM)
-gmx_test_inline_asm_gcc_x86(GMX_X86_GCC_INLINE_ASM)
-gmx_test_inline_asm_msvc_x86(GMX_X86_MSVC_INLINE_ASM)
-
-# turn on SSE if supported with reasonable defaults.
-if (${GMX_ACCELERATION} STREQUAL "auto" AND NOT GMX_OPENMM)
- if(CMAKE_SYSTEM_PROCESSOR MATCHES "^(i.86|x86|x64|x86_64|AMD64|amd64)" OR CYGWIN)
-
- set(GMX_ACCELERATION "SSE" CACHE STRING "Accelerated kernels. Pick one of: auto, none, SSE, BlueGene, Power6, ia64, altivec, fortran" FORCE)
-
- # Determine the assembler/compiler to use
- else()
- set(GMX_ACCELERATION "none" CACHE STRING "Accelerated kernels. Pick one of: auto, none, SSE, BlueGene, Power6, ia64, altivec, fortran" FORCE)
- endif()
-endif (${GMX_ACCELERATION} STREQUAL "auto" AND NOT GMX_OPENMM)
-
include(gmxTestXDR)
gmx_test_xdr(GMX_SYSTEM_XDR)
if(NOT GMX_SYSTEM_XDR)
# Process nonbonded accelerated kernels settings
string(TOUPPER ${GMX_ACCELERATION} ${GMX_ACCELERATION})
if(${GMX_ACCELERATION} STREQUAL "NONE")
- # nothing to do
-elseif(${GMX_ACCELERATION} STREQUAL "SSE")
-
- if (GMX_64_BIT)
- set(GMX_X86_64_ASM ON CACHE BOOL "Add SSE assembly files for x86_64" FORCE)
- else (GMX_64_BIT)
- set(GMX_IA32_ASM ON CACHE BOOL "Add SSE assembly files for i386" FORCE)
- endif (GMX_64_BIT)
-
- if( WIN32 AND NOT CYGWIN )
- option(GMX_ASM_USEASM_NASM "Use Nasm for assembly instead of compiler (necessary on windows)" ON)
- else()
- option(GMX_ASM_USEASM_NASM "Use Nasm for assembly instead of compiler (necessary on windows)" OFF)
+ # nothing to do
+elseif(${GMX_ACCELERATION} STREQUAL "SSE2")
+
+ GMX_TEST_CFLAG(GNU_SSE2_CFLAG "-msse2" GROMACS_C_FLAGS)
+ if(NOT GNU_SSE2_CFLAG)
+ GMX_TEST_CFLAG(MSVC_SSE2_CFLAG "/arch:SSE2" GROMACS_C_FLAGS)
+ endif(NOT GNU_SSE2_CFLAG)
+
+ GMX_TEST_CXXFLAG(GNU_SSE2_CXXFLAG "-msse2" GROMACS_CXX_FLAGS)
+ if(NOT GNU_SSE2_CXXFLAG)
+ GMX_TEST_CXXFLAG(MSVC_SSE2_CXXFLAG "/arch:SSE2" GROMACS_CXX_FLAGS)
+ endif(NOT GNU_SSE2_CXXFLAG)
+
+ # We dont warn for lacking SSE2 flag support, since that is probably standard today.
+
+ # Only test the include after we have tried to add the correct flag for SSE2 support
+ check_include_file(emmintrin.h HAVE_EMMINTRIN_H ${GROMACS_C_FLAGS})
+
+ if(NOT HAVE_EMMINTRIN_H)
+ message(FATAL_ERROR "Cannot find emmintrin.h, which is required for SSE2 intrinsics support.")
+ endif(NOT HAVE_EMMINTRIN_H)
+
+ # The user should not be able to set this orthogonally to the acceleration
+ set(GMX_X86_SSE2 1)
+ if (NOT ACCELERATION_QUIETLY)
+ message(STATUS "Enabling SSE2 Gromacs acceleration, and it will help compiler optimization.")
endif()
- mark_as_advanced(GMX_ASM_USEASM_NASM)
-
- if (NOT GMX_64_BIT)
- # for 32-bit compiles, we might need to turn on sse
- CHECK_C_COMPILER_FLAG("-msse2" XFLAGS_SSE)
- if (XFLAGS_SSE)
- set(GROMACS_C_FLAGS "-msse2 ${GROMACS_C_FLAGS}")
- endif (XFLAGS_SSE)
- CHECK_CXX_COMPILER_FLAG("-msse2" XXFLAGS_SSE)
- if (XXFLAGS_SSE)
- set(GROMACS_CXX_FLAGS "-msse2 ${GROMACS_CXX_FLAGS}")
- endif (XXFLAGS_SSE)
- endif (NOT GMX_64_BIT)
-
- if(GMX_IA32_ASM)
- if(GMX_DOUBLE)
- set(GMX_IA32_SSE2 1)
- else()
- set(GMX_IA32_SSE 1)
- endif()
- elseif(GMX_X86_64_ASM)
- if(GMX_DOUBLE)
- set(GMX_X86_64_SSE2 1)
- else()
- set(GMX_X86_64_SSE 1)
- endif()
+
+elseif(${GMX_ACCELERATION} STREQUAL "SSE4.1")
+
+ GMX_TEST_CFLAG(GNU_SSE4_CFLAG "-msse4.1" GROMACS_C_FLAGS)
+ if (NOT GNU_SSE4_CFLAG)
+ GMX_TEST_CFLAG(MSVC_SSE4_CFLAG "/arch:SSE4.1" GROMACS_C_FLAGS)
+ endif(NOT GNU_SSE4_CFLAG)
+ if (NOT GNU_SSE4_CFLAG AND NOT MSVC_SSE4_CFLAG)
+ message(WARNING "No C SSE4.1 flag found. Consider a newer compiler, or disable SSE4.1 for slightly lower performance.")
+ # Not surprising if we end up here! MSVC current does not support the SSE4.1 flag. However, it appears to accept SSE4.1
+ # intrinsics when SSE2 support is enabled, so we try that instead.
+ GMX_TEST_CFLAG(MSVC_SSE2_CFLAG "/arch:SSE2" GROMACS_C_FLAGS)
+ endif(NOT GNU_SSE4_CFLAG AND NOT MSVC_SSE4_CFLAG)
+
+ GMX_TEST_CXXFLAG(GNU_SSE4_CXXFLAG "-msse4.1" GROMACS_CXX_FLAG)
+ if (NOT GNU_SSE4_CXXFLAG)
+ GMX_TEST_CXXFLAG(MSVC_SSE4_CXXFLAG "/arch:SSE4.1" GROMACS_CXX_FLAGS)
+ endif(NOT GNU_SSE4_CXXFLAG)
+ if (NOT GNU_SSE4_CXXFLAG AND NOT MSVC_SSE4_CXXFLAG)
+ message(WARNING "No C++ SSE4.1 flag found. Consider a newer compiler, or disable SSE4.1 for slightly lower performance.")
+ # Not surprising if we end up here! MSVC current does not support the SSE4.1 flag. However, it appears to accept SSE4.1
+ # intrinsics when SSE2 support is enabled, so we try that instead.
+ GMX_TEST_CXXFLAG(MSVC_SSE2_CXXFLAG "/arch:SSE2" GROMACS_CXX_FLAGS)
+ endif(NOT GNU_SSE4_CXXFLAG AND NOT MSVC_SSE4_CXXFLAG)
+
+ # This must come after we have added the -msse4.1 flag on some platforms.
+ check_include_file(smmintrin.h HAVE_SMMINTRIN_H ${GROMACS_C_FLAGS})
+
+ if(NOT HAVE_SMMINTRIN_H)
+ message(FATAL_ERROR "Cannot find smmintrin.h, which is required for SSE4.1 intrinsics support.")
+ endif(NOT HAVE_SMMINTRIN_H)
+
+ # The user should not be able to set this orthogonally to the acceleration
+ set(GMX_X86_SSE4_1 1)
+ set(GMX_X86_SSE2 1)
+ if (NOT ACCELERATION_QUIETLY)
+ message(STATUS "Enabling SSE4.1 Gromacs acceleration, and it will help compiler optimization.")
endif()
-# if(HAVE_PMMINTRIN_H)
-# set(GMX_SSE3 1)
-# endif(HAVE_PMMINTRIN_H)
-# if(HAVE_SMMINTRIN_H)
-# set(GMX_SSE4_1 1)
-# endif(HAVE_SMMINTRIN_H)
+elseif(${GMX_ACCELERATION} STREQUAL "AVX_128_FMA" OR ${GMX_ACCELERATION} STREQUAL "AVX_256")
+
+ # Set the AVX compiler flag for both these choices!
+
+ GMX_TEST_CFLAG(GNU_AVX_CFLAG "-mavx" GROMACS_C_FLAGS)
+ if (NOT GNU_AVX_CFLAG)
+ GMX_TEST_CFLAG(MSVC_AVX_CFLAG "/arch:AVX" GROMACS_C_FLAGS)
+ endif (NOT GNU_AVX_CFLAG)
+ if (NOT GNU_AVX_CFLAG AND NOT MSVC_AVX_CFLAG)
+ message(WARNING "No C AVX flag found. Consider a newer compiler, or disable AVX for much lower performance.")
+ endif (NOT GNU_AVX_CFLAG AND NOT MSVC_AVX_CFLAG)
+
+ GMX_TEST_CXXFLAG(GNU_AVX_CXXFLAG "-mavx" GROMACS_CXX_FLAGS)
+ if (NOT GNU_AVX_CXXFLAG)
+ GMX_TEST_CXXFLAG(MSVC_AVX_CXXFLAG "/arch:AVX" GROMACS_CXX_FLAGS)
+ endif (NOT GNU_AVX_CXXFLAG)
+ if (NOT GNU_AVX_CXXFLAG AND NOT MSVC_AVX_CXXFLAG)
+ message(WARNING "No C++ AVX flag found. Consider a newer compiler, or disable AVX for much lower performance.")
+ endif (NOT GNU_AVX_CXXFLAG AND NOT MSVC_AVX_CXXFLAG)
+
+ # Only test the header after we have tried to add the flag for AVX support
+ check_include_file(immintrin.h HAVE_IMMINTRIN_H ${GROMACS_C_FLAGS})
+
+ if(NOT HAVE_IMMINTRIN_H)
+ message(FATAL_ERROR "Cannot find immintrin.h, which is required for AVX intrinsics support. Consider switching compiler.")
+ endif(NOT HAVE_IMMINTRIN_H)
+
+ # AMD says we should include x86intrin.h for FMA support, but MSVC seems to do fine without it, so don't require it.
+ check_include_file(x86intrin.h HAVE_X86INTRIN_H ${GROMACS_C_FLAGS})
+
+ # The user should not be able to set this orthogonally to the acceleration
+ set(GMX_X86_SSE4_1 1)
+ set(GMX_X86_SSE2 1)
+
+ # But just enable one of the choices internally...
+ if(${GMX_ACCELERATION} STREQUAL "AVX_128_FMA")
+ set(GMX_X86_AVX_128_FMA 1)
+ if (NOT ACCELERATION_QUIETLY)
+ message(STATUS "Enabling 128-bit AVX Gromacs acceleration (with fused-multiply add), and it will help compiler optimization.")
+ endif()
+ else()
+ # If we are not doing AVX_128, it must be AVX_256...
+ set(GMX_X86_AVX_256 1)
+ if (NOT ACCELERATION_QUIETLY)
+ message(STATUS "Enabling 256-bit AVX Gromacs acceleration, and it will help compiler optimization.")
+ endif()
+ endif()
elseif(${GMX_ACCELERATION} STREQUAL "FORTRAN")
- set(GMX_FORTRAN 1)
- #these are switch on by default sometimes
- set(GMX_IA32_ASM 0)
- set(GMX_GMX_X86_64_ASM 0)
+
+# Fortran is temporarily disabled while we push in nbNxN kernels.
+# We need to fake it a bit here to avoid jenkins build errors!
+# add_definitions(-DGMX_FORTRAN)
+
elseif(${GMX_ACCELERATION} STREQUAL "BLUEGENE")
# GMX_ACCELERATION=BlueGene should be set in the Toolchain-BlueGene?-???.cmake file
- message(STATUS "Configuring for BlueGene")
+ if (NOT ACCELERATION_QUIETLY)
+ message(STATUS "Configuring for BlueGene")
+ endif()
set(GMX_BLUEGENE 1)
if (${CMAKE_SYSTEM_NAME} STREQUAL "BlueGeneL")
set(SHARED_LIBS_DEFAULT OFF CACHE BOOL "Shared libraries not compatible with BlueGene/L, disabled!" FORCE)
set(GMX_SOFTWARE_INVSQRT OFF CACHE BOOL "Do not use software reciprocal square root on BlueGene" FORCE)
set(GMX_POWERPC_INVSQRT ON CACHE BOOL "Use hardware reciprocal square root on BlueGene" FORCE)
set(GMX_X11 OFF CACHE BOOL "X11 not compatible with BlueGene, disabled!" FORCE)
- set(GMX_THREADS OFF CACHE BOOL "Threads not compatible with BlueGene, disabled!" FORCE)
+ set(GMX_THREAD_MPI OFF CACHE BOOL "Thread-MPI not compatible with BlueGene, disabled!" FORCE)
set(GMX_MPI ON CACHE BOOL "Use MPI on BlueGene" FORCE)
+ # Access to /etc/passwd is not available on the back end of BlueGene,
+ # despite being detected by CMake. This can cause linker warnings
+ # about harmless things in src/gmxlib/string2.h.
+ set(HAVE_PWD_H OFF)
+ # The automatic testing for endianness does not work for the BlueGene cross-compiler
+ set(GMX_IEEE754_BIG_ENDIAN_BYTE_ORDER 1 CACHE INTERNAL "BlueGene has big endian FP byte order (by default)" FORCE)
+ set(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER 1 CACHE INTERNAL "BlueGene has big endian FP word order (by default)" FORCE)
elseif(${GMX_ACCELERATION} STREQUAL "POWER6")
set(GMX_POWER6 1)
set(GMX_SOFTWARE_INVSQRT OFF CACHE BOOL "Do not use software reciprocal square root on Power6" FORCE)
set(GMX_POWERPC_INVSQRT ON CACHE BOOL "Use hardware reciprocal square root on Power6" FORCE)
-elseif(${GMX_ACCELERATION} STREQUAL "IA64")
- set(GMX_IA64_ASM 1)
- set(DISABLE_WATERWATER_NLIST 1)
- set(DISABLE_WATER_NLIST 1)
-elseif(${GMX_ACCELERATION} STREQUAL "ALTIVEC")
- check_include_files(altivec.h HAVE_ALTIVEC_H)
- if(HAVE_ALTIVEC_H)
- set(GMX_PPC_ALTIVEC 1)
- endif(HAVE_ALTIVEC_H)
else(${GMX_ACCELERATION} STREQUAL "NONE")
- MESSAGE(FATAL_ERROR "Unrecognized option for accelerated kernels: ${GMX_ACCELERATION}. Pick one of auto, none, SSE, Fortran, BlueGene, Power6, ia64, altivec")
+ MESSAGE(FATAL_ERROR "Unrecognized option for accelerated kernels: ${GMX_ACCELERATION}. Pick one of None, SSE2, SSE4.1, AVX_128_FMA, AVX_256, Fortran, BlueGene, Power6")
endif(${GMX_ACCELERATION} STREQUAL "NONE")
+set(ACCELERATION_QUIETLY TRUE CACHE INTERNAL "")
if(GMX_FORTRAN OR GMX_POWER6)
- if (GMX_THREADS)
- message(FATAL_ERROR "FORTRAN/POWER6 is incompatible with threads and only provides a speed-up on certain IBM compilers. Disable FORTRAN (or threads if you really want to use FORTRAN kernels).")
- endif(GMX_THREADS)
+ if (GMX_THREAD_MPI)
+ message(FATAL_ERROR "FORTRAN/POWER6 is incompatible with thread-MPI and only provides a speed-up on certain IBM compilers. Disable FORTRAN (or threads if you really want to use FORTRAN kernels).")
+ endif(GMX_THREAD_MPI)
enable_language(Fortran)
include(FortranCInterface)
discover_fortran_mangling(prefix isupper suffix extra_under_score found)
set(PKG_FFT "")
set(PKG_FFT_LIBS "")
if(${GMX_FFT_LIBRARY} STREQUAL "FFTW3")
-# MESSAGE(STATUS "Using external FFT library - fftw3")
if(GMX_DOUBLE)
- find_package(FFTW3 REQUIRED)
- include_directories(${FFTW3_INCLUDE_DIR})
- set(FFT_LIBRARIES ${FFTW3_LIBRARIES})
- set(PKG_FFT "fftw3")
+ find_package(FFTW 3 COMPONENTS fftw)
else(GMX_DOUBLE)
- find_package(FFTW3F REQUIRED)
- include_directories(${FFTW3F_INCLUDE_DIR})
- set(FFT_LIBRARIES ${FFTW3F_LIBRARIES})
- set(PKG_FFT "fftw3f")
+ find_package(FFTW 3 COMPONENTS fftwf)
endif(GMX_DOUBLE)
-
- if(NOT FFTW3_FOUND AND NOT FFTW3F_FOUND)
- MESSAGE(FATAL_ERROR "Cannot find fftw3 (with correct precision). Fix it, choose another FFT library, or use the Gromacs built-in fftpack (slower)!")
- endif(NOT FFTW3_FOUND AND NOT FFTW3F_FOUND)
+ if(NOT FFTW_FOUND)
+ MESSAGE(FATAL_ERROR "Cannot find FFTW3 (with correct precision - libfftw3f for single precision GROMACS or libfftw3 for double precision GROMACS). Fix it, choose another FFT library, or use the Gromacs built-in fftpack (slower)!")
+ endif(NOT FFTW_FOUND)
+ include_directories(${FFTW_INCLUDE_DIRS})
+ set(FFT_LIBRARIES ${FFTW_LIBRARIES})
+ set(PKG_FFT "${FFTW_PKG}")
set(GMX_FFT_FFTW3 1)
-elseif(${GMX_FFT_LIBRARY} STREQUAL "FFTW2")
-# MESSAGE(STATUS "Using external FFT library - fftw2")
- if(GMX_DOUBLE)
- find_package(FFTW2 REQUIRED)
- else(GMX_DOUBLE)
- find_package(FFTW2F REQUIRED)
- endif(GMX_DOUBLE)
+ if (${GMX_ACCELERATION} STREQUAL "SSE" AND NOT FFTW_HAVE_SSE)
+ message(WARNING "The fftw library found is compiled without SSE support, which makes it slow. Consider recompiling it or contact your admin")
+ endif (${GMX_ACCELERATION} STREQUAL "SSE" AND NOT FFTW_HAVE_SSE)
- if(NOT FFTW2_FOUND)
- MESSAGE(FATAL_ERROR "Cannot find fftw2 (with correct precision). Fix it, choose another FFT library, or use the Gromacs built-in fftpack (slower)!")
- endif(NOT FFTW2_FOUND)
- include_directories(${FFTW2_INCLUDE_DIR})
- set(FFT_LIBRARIES ${FFTW2_LIBRARIES})
- set(PKG_FFT_LIBS ${FFTW2_LIBRARIES})
-
- if("${FFTW2_LIBRARIES}" MATCHES "dfftw")
- set(FFTW2_NAME_DFFTW 1)
- elseif("${FFTW2_LIBRARIES}" MATCHES "sfftw")
- set(FFTW2_NAME_SFFTW 1)
- else("${FFTW2_LIBRARIES}" MATCHES "dfftw")
- set(FFTW2_NAME_FFTW 1)
- endif("${FFTW2_LIBRARIES}" MATCHES "dfftw")
-
- set(GMX_FFT_FFTW2 1)
elseif(${GMX_FFT_LIBRARY} STREQUAL "MKL")
# MESSAGE(STATUS "Using external FFT library - Intel MKL")
find_package(MKL REQUIRED)
MESSAGE(STATUS "Using internal FFT library - fftpack")
set(GMX_FFT_FFTPACK 1)
else(${GMX_FFT_LIBRARY} STREQUAL "FFTW3")
- MESSAGE(FATAL_ERROR "Invalid FFT library setting: ${GMX_FFT_LIBRARY}. Choose one of: fftw3, fftw2, mkl, acml, fftpack")
+ MESSAGE(FATAL_ERROR "Invalid FFT library setting: ${GMX_FFT_LIBRARY}. Choose one of: fftw3, mkl, fftpack")
endif(${GMX_FFT_LIBRARY} STREQUAL "FFTW3")
+set(GMX_EXTERNAL_BLAS TRUE CACHE BOOL "Use external BLAS instead of built-in")
+set(GMX_EXTERNAL_LAPACK TRUE CACHE BOOL "Use external LAPACK instead of built-in")
# MKL has BLAS/LAPACK routines
-if(HAVE_MKL OR ACCELERATE_FRAMEWORK)
- set(GMX_EXTERNAL_BLAS TRUE CACHE BOOL "Use external BLAS instead of built-in")
- set(GMX_EXTERNAL_LAPACK TRUE CACHE BOOL "Use external LAPACK instead of built-in")
-else(HAVE_MKL OR ACCELERATE_FRAMEWORK)
- set(GMX_EXTERNAL_BLAS FALSE CACHE BOOL "Use external BLAS instead of built-in")
- set(GMX_EXTERNAL_LAPACK FALSE CACHE BOOL "Use external LAPACK instead of built-in")
+if(NOT HAVE_MKL AND NOT ACCELERATE_FRAMEWORK)
if(GMX_EXTERNAL_BLAS)
- #FindBLAS needs Fortran
- enable_language(Fortran)
if (GMX_BLAS_USER)
list(APPEND GMX_EXTRA_LIBRARIES ${GMX_BLAS_USER})
else(GMX_BLAS_USER)
- find_package(BLAS REQUIRED)
- list(APPEND GMX_EXTRA_LIBRARIES ${BLAS_LIBRARIES})
+ set(BLAS_FIND_QUIETLY ON)
+ find_package(BLAS)
+ if (BLAS_FOUND)
+ list(APPEND GMX_EXTRA_LIBRARIES ${BLAS_LIBRARIES})
+ else()
+ MESSAGE(STATUS "Using internal BLAS library")
+ set(GMX_EXTERNAL_BLAS FALSE CACHE BOOL "Use external BLAS instead of built-in" FORCE)
+ endif()
endif(GMX_BLAS_USER)
endif(GMX_EXTERNAL_BLAS)
if(GMX_EXTERNAL_LAPACK)
- #FindLAPACK needs Fortran
- enable_language(Fortran)
- if (GMX_LAPACK_USER)
+ if (GMX_LAPACK_USER)
list(APPEND GMX_EXTRA_LIBRARIES ${GMX_LAPACK_USER})
- else(GMX_LAPACK_USER)
- find_package(LAPACK REQUIRED)
- list(APPEND GMX_EXTRA_LIBRARIES ${LAPACK_LIBRARIES})
+ else(GMX_LAPACK_USER)
+ set(LAPACK_FIND_QUIETLY ON)
+ find_package(LAPACK)
+ if (LAPACK_FOUND)
+ list(APPEND GMX_EXTRA_LIBRARIES ${LAPACK_LIBRARIES})
+ else()
+ MESSAGE(STATUS "Using internal LAPACK library")
+ set(GMX_EXTERNAL_LAPACK FALSE CACHE BOOL "Use external LAPACK instead of built-in" FORCE)
+ endif()
endif(GMX_LAPACK_USER)
endif(GMX_EXTERNAL_LAPACK)
-endif(HAVE_MKL OR ACCELERATE_FRAMEWORK)
+endif()
mark_as_advanced(GMX_EXTERNAL_LAPACK)
mark_as_advanced(GMX_EXTERNAL_BLAS)
+set(GMX_USE_PLUGINS OFF CACHE INTERNAL "Whether GROMACS will really try to compile support for VMD plugins")
+set(GMX_VMD_PLUGIN_PATH)
+mark_as_advanced(GMX_VMD_PLUGIN_PATH)
+
+if(GMX_LOAD_PLUGINS)
+ if(CYGWIN OR NOT WIN32)
+ # Native Windows does not have, nor need dlopen
+ # Note that WIN32 is set with Cygwin, but Cygwin needs dlopen to use plug-ins
+ include(gmxTestdlopen)
+ gmx_test_dlopen(HAVE_DLOPEN)
+ endif()
+
+ find_package(VMD)
+
+ # Test for unsuitable versions of VMD
+
+ if(VMD_FOUND AND NOT GMX_VMD_PLUGIN_PATH)
+ message(STATUS "Checking for suitable VMD version")
+ exec_program(${VMD_EXECUTABLE}
+ ARGS --help
+ OUTPUT_VARIABLE VMD_HELP
+ RETURN_VALUE VMD_EXEC_RETURN)
+
+ if(VMD_EXEC_RETURN EQUAL 0)
+ # This is the accepted idiom for subexpression matching, unfortunately
+ string(REGEX REPLACE ".*VMD for .*, version ([0-9]+\\.[0-9]*\\.?[0-9]*).*" "\\1" VMD_VERSION ${VMD_HELP})
+ string(REGEX REPLACE ".*VMD for (.*), version .*" "\\1" VMD_ARCH ${VMD_HELP})
+
+ if(VMD_VERSION VERSION_LESS "1.8")
+ MESSAGE(WARNING "Found VMD version ${VMD_VERSION}, but GROMACS needs at least 1.8")
+ unset(VMD_EXECUTABLE)
+ set(VMD_FOUND FALSE)
+ else()
+ message(STATUS "VMD version ${VMD_VERSION} is suitable")
+ if(DEFINED ENV{VMDDIR})
+ # This permits GROMACS to avoid hard-coding a fall-back
+ # path that it can tell right now would be useless.
+ set(GMX_VMD_PLUGIN_PATH "$ENV{VMDDIR}/plugins/${VMD_ARCH}/molfile" CACHE PATH "Path to VMD plugins for molfile I/O")
+ else()
+ set(GMX_VMD_PLUGIN_PATH "/usr/local/lib/vmd/plugins/*/molfile" CACHE PATH "Path to VMD plugins for molfile I/O")
+ endif()
+ endif()
+
+ # clean up
+ unset(VMD_HELP)
+ unset(VMD_VERSION)
+ unset(VMD_ARCH)
+ endif()
+ endif()
+
+ # so, should we use plug-ins?
+ if((WIN32 AND NOT CYGWIN) OR (HAVE_DLOPEN AND BUILD_SHARED_LIBS))
+ if (NOT VMD_QUIETLY)
+ MESSAGE(STATUS "Found the ability to use plug-ins when building shared libaries, so will compile to use plug-ins (e.g. to read VMD-supported file formats).")
+ endif(NOT VMD_QUIETLY)
+ set(GMX_USE_PLUGINS ON)
+ list(APPEND GMX_EXTRA_LIBRARIES ${CMAKE_DL_LIBS}) # magic cross-platform pre-set variable for dlopen library
+ set(PKG_DL_LIBS "-l${CMAKE_DL_LIBS}")
+ else()
+ set(PKG_DL_LIBS)
+ endif()
+endif(GMX_LOAD_PLUGINS)
+set(VMD_QUIETLY TRUE CACHE INTERNAL "")
+
+if(EXISTS "${CMAKE_SOURCE_DIR}/admin/.isreposource")
+ if(NOT CMAKE_CROSSCOMPILING)
+ option(GMX_BUILD_MANPAGES "Build man pages" ON)
+ else()
+ message(STATUS "Building the man pages is not available when "
+ "cross-compiling the developer version from git")
+ endif()
+else()
+ #make sure source package contains all man pages
+ if(NOT EXISTS "${CMAKE_SOURCE_DIR}/man/man1/ngmx.1")
+ message(FATAL_ERROR "Man pages are missing from source package.")
+ endif()
+endif()
+mark_as_advanced(GMX_BUILD_MANPAGES)
+
# Math and thread libraries must often come after all others when linking...
if(HAVE_LIBM)
list(APPEND GMX_EXTRA_LIBRARIES m)
CACHE STRING "Linker flags" FORCE)
endif (NOT DEFINED GROMACS_C_FLAGS_SET)
+######################################
+# Output compiler and CFLAGS used
+######################################
+execute_process(COMMAND ${CMAKE_C_COMPILER} --version RESULT_VARIABLE TMP_RESULT OUTPUT_VARIABLE CC_VERSION ERROR_VARIABLE CC_VERSION)
+#try executing just the compiler command if that failed...
+if(TMP_RESULT)
+ execute_process(COMMAND ${CMAKE_C_COMPILER} RESULT_VARIABLE TMP_RESULT OUTPUT_VARIABLE CC_VERSION ERROR_VARIABLE CC_VERSION)
+endif()
+string(LENGTH "${CC_VERSION}" len)
+if(len)
+ string(REGEX MATCH "[^\n]*" CC_VERSION "${CC_VERSION}")
+endif()
+
+set(BUILD_COMPILER "@CMAKE_C_COMPILER@ @CMAKE_C_COMPILER_ID@ @CC_VERSION@")
+set(BUILD_CFLAGS "@CMAKE_C_FLAGS@")
+if(CMAKE_BUILD_TYPE STREQUAL "Debug")
+ set(BUILD_CFLAGS "@BUILD_CFLAGS@ @CMAKE_C_FLAGS_DEBUG@")
+elseif(CMAKE_BUILD_TYPE STREQUAL "Release")
+ set(BUILD_CFLAGS "@BUILD_CFLAGS@ @CMAKE_C_FLAGS_RELEASE@")
+elseif(CMAKE_BUILD_TYPE STREQUAL "MinSizeRel")
+ set(BUILD_CFLAGS "@BUILD_CFLAGS@ @CMAKE_C_FLAGS_MINSIZEREL@")
+elseif(CMAKE_BUILD_TYPE STREQUAL "RelWithDebInfo")
+ set(BUILD_CFLAGS "@BUILD_CFLAGS@ @CMAKE_C_FLAGS_RELWITHDEBINFO@")
+endif()
+
+
########################################################################
# Specify install locations and which subdirectories to process #
########################################################################
-if ( DEFINED GMXLIB )
- set(LIB_INSTALL_DIR "${CMAKE_INSTALL_PREFIX}/${GMXLIB}")
+if (GMX_USE_RELATIVE_INSTALL_PATH)
+ set(GMX_INSTALL_PREFIX "" CACHE STRING "Prefix gets appended to CMAKE_INSTALL_PREFIX. For cpack it sets the root folder of the archive.")
+ mark_as_advanced(GMX_INSTALL_PREFIX)
else()
- set(LIB_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/lib)
+ set(GMX_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}/")
+endif()
+
+if ( NOT DEFINED GMXLIB )
+ set(GMXLIB lib)
endif()
-set(BIN_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/bin)
-set(DATA_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/share/gromacs)
-set(MAN_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/share/man)
-set(INCL_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/include)
+set(LIB_INSTALL_DIR "${GMX_INSTALL_PREFIX}${GMXLIB}")
+set(BIN_INSTALL_DIR ${GMX_INSTALL_PREFIX}bin)
+set(DATA_INSTALL_DIR ${GMX_INSTALL_PREFIX}share/gromacs)
+set(MAN_INSTALL_DIR ${GMX_INSTALL_PREFIX}share/man)
+set(INCL_INSTALL_DIR ${GMX_INSTALL_PREFIX}include)
set(GMXLIBDIR ${DATA_INSTALL_DIR}/top)
if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
set(CMAKE_SKIP_BUILD_RPATH FALSE)
set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
- set(CMAKE_INSTALL_RPATH "${LIB_INSTALL_DIR}")
+ set(CMAKE_INSTALL_RPATH "\\\$ORIGIN/../${GMXLIB}")
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
endif()
#define int64_t long long
#endif
-#if (defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64 || defined __CYGWIN__ || defined __CYGWIN32__ || defined GMX_INTERNAL_XDR)
+#include "gmx_header_config.h"
+#if (defined GMX_NATIVE_WINDOWS || defined GMX_CYGWIN || defined GMX_INTERNAL_XDR)
#include "gmx_system_xdr.h"
#else
#include <rpc/rpc.h>
float xdr_xtc_estimate_dt(FILE *fp, XDR *xdrs, int natoms, gmx_bool * bOK);
- int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms);
+ int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms,gmx_bool bSeekForwardOnly);
int xdr_xtc_seek_frame(int frame, FILE *fp, XDR *xdrs, int natoms);
#include "gmxfio.h"
#include "md5.h"
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "thread_mpi.h"
#endif
static t_fileio *open_files = NULL;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* this mutex locks the open_files structure so that no two threads can
modify it.
type of access to the fio's elements. */
void gmx_fio_lock(t_fileio *fio)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Lock_lock(&(fio->mtx));
#endif
}
/* unlock the mutex associated with this fio. */
void gmx_fio_unlock(t_fileio *fio)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Lock_unlock(&(fio->mtx));
#endif
}
open_files->fn=NULL;
open_files->next=open_files;
open_files->prev=open_files;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Lock_init(&(open_files->mtx));
#endif
}
static void gmx_fio_insert(t_fileio *fio)
{
t_fileio *prev;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* first lock the big open_files mutex. */
tMPI_Thread_mutex_lock(&open_file_mutex);
#endif
gmx_fio_unlock(open_files);
gmx_fio_unlock(fio);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* now unlock the big open_files mutex. */
tMPI_Thread_mutex_unlock(&open_file_mutex);
#endif
t_fileio *ret;
/* first lock the big open_files mutex and the dummy's mutex */
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* first lock the big open_files mutex. */
tMPI_Thread_mutex_lock(&open_file_mutex);
#endif
if (fio->next==open_files)
{
ret=NULL;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&open_file_mutex);
#endif
}
static void gmx_fio_stop_getting_next(t_fileio *fio)
{
gmx_fio_unlock(fio);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&open_file_mutex);
#endif
}
}
snew(fio, 1);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Lock_init(&(fio->mtx));
#endif
bRead = (newmode[0]=='r' && newmode[1]!='+');
{
int rc = 0;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* first lock the big open_files mutex. */
/* We don't want two processes operating on the list at the same time */
tMPI_Thread_mutex_lock(&open_file_mutex);
sfree(fio);
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&open_file_mutex);
#endif
return ret;
}
- int xtc_seek_time(t_fileio *fio, real time, int natoms)
+ int xtc_seek_time(t_fileio *fio, real time, int natoms,gmx_bool bSeekForwardOnly)
{
int ret;
gmx_fio_lock(fio);
- ret=xdr_xtc_seek_time(time, fio->fp, fio->xdr, natoms);
+ ret=xdr_xtc_seek_time(time, fio->fp, fio->xdr, natoms, bSeekForwardOnly);
gmx_fio_unlock(fio);
return ret;
static XDR *xdridptr[MAXID];
static char xdrmodes[MAXID];
static unsigned int cnt;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* we need this because of the global variables above for FORTRAN binding.
The I/O operations are going to be slow. */
static tMPI_Thread_mutex_t xdr_fortran_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
static void xdr_fortran_lock(void)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&xdr_fortran_mutex);
#endif
}
static void xdr_fortran_unlock(void)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&xdr_fortran_mutex);
#endif
}
char newtype[5];
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
if (!tMPI_Thread_mutex_trylock( &xdr_fortran_mutex ))
{
tMPI_Thread_mutex_unlock( &xdr_fortran_mutex );
int xdrid;
int rc = 0;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
if (!tMPI_Thread_mutex_trylock( &xdr_fortran_mutex ))
{
tMPI_Thread_mutex_unlock( &xdr_fortran_mutex );
int bytes[32];
int i, j, num_of_bytes, p, num;
- bytes[1] = bytes[2] = bytes[3] = 0;
+ bytes[0] = bytes[1] = bytes[2] = bytes[3] = 0;
num_of_bytes = 0;
while (num_of_bits > 8) {
bytes[num_of_bytes++] = receivebits(buf, 8);
- int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms)
+ int xdr_xtc_seek_time(real time, FILE *fp, XDR *xdrs, int natoms,gmx_bool bSeekForwardOnly)
{
float t;
float dt;
int res;
int dt_sign = 0;
+ if (bSeekForwardOnly)
+ {
+ low = gmx_ftell(fp);
+ }
if (gmx_fseek(fp,0,SEEK_END))
{
return -1;
/* round to int */
high /= XDR_INT_SIZE;
high *= XDR_INT_SIZE;
- offset = ((high / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
+ offset = (((high-low) / 2) / XDR_INT_SIZE) * XDR_INT_SIZE;
if (gmx_fseek(fp,offset,SEEK_SET))
{
* ====================================================
*/
-#if (INT_MAX == 2147483647)
+#if ( (defined SIZEOF_INT && SIZEOF_INT==4) || (SIZEOF_INT_MAX == 2147483647) )
typedef int erf_int32_t;
typedef unsigned int erf_u_int32_t;
#elif (LONG_MAX == 2147483647L)
conv.d=x;
- #ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
+ /* In release-4-6 and later branches, only the test for
+ * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
+ #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
hx=conv.i[0];
#else
hx=conv.i[1];
conv.d = x;
- #ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
+ /* In release-4-6 and later branches, only the test for
+ * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
+ #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
conv.i[1] = 0;
#else
conv.i[0] = 0;
conv.d = x;
- #ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
+ /* In release-4-6 and later branches, only the test for
+ * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
+ #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
hx=conv.i[0];
#else
hx=conv.i[1];
conv.d = x;
- #ifdef IEEE754_BIG_ENDIAN_WORD_ORDER
+ /* In release-4-6 and later branches, only the test for
+ * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
+ #if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
conv.i[1] = 0;
#else
conv.i[0] = 0;
#endif
-float fast_float_erf(float x)
-{
- float t,ans;
-
- t=1.0/(1.0+0.5*x);
- ans=t*exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
- t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
- t*(-0.82215223+t*0.17087277)))))))));
- return 1.0-ans;
-}
-
-float fast_float_erfc(float x)
-{
- float t,ans;
-
- t=1.0/(1.0+0.5*x);
- ans=t*exp(-x*x-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
- t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
- t*(-0.82215223+t*0.17087277)))))))));
- return ans;
-}
-
gmx_bool gmx_isfinite(real x)
{
gmx_bool returnval = TRUE;
#include <smalloc.h>
#include <string2.h>
#include <vec.h>
+#include <assert.h>
#include <indexutil.h>
#include <poscalc.h>
static void
init_item_evaloutput(t_selelem *sel)
{
+ assert(!(sel->child == NULL &&
+ (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)));
+
/* Process children. */
if (sel->type != SEL_SUBEXPRREF)
{
&& ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
|| (sel->cdata->flags & SEL_CDATA_FULLEVAL)))
{
+ assert(sel->child);
sel->cdata->gmin = sel->child->cdata->gmin;
sel->cdata->gmax = sel->child->cdata->gmax;
}
{
gmx_ana_index_reserve(sel->child->v.u.g, gmin->isize);
gmx_ana_index_copy(sel->child->v.u.g, gmin, FALSE);
- if (sel->child->u.cgrp.isize > 0)
+ if (sel->child->u.cgrp.nalloc_index > 0)
{
+ /* Keep the name as in evaluate_boolean_static_part(). */
+ char *name = sel->child->u.cgrp.name;
gmx_ana_index_reserve(&sel->child->u.cgrp, gmin->isize);
gmx_ana_index_copy(&sel->child->u.cgrp, gmin, FALSE);
+ sel->child->u.cgrp.name = name;
+ }
+ else
+ {
+ sel->child->u.cgrp.isize = sel->child->v.u.g->isize;
}
}
break;
case SEL_EXPRESSION:
case SEL_MODIFIER:
+ assert(g);
rc = _gmx_sel_evaluate_method_params(data, sel, g);
if (rc != 0)
{
}
else if (sel->u.cgrp.isize == 0)
{
+ assert(g);
gmx_ana_index_reserve(&sel->u.cgrp, g->isize);
rc = sel->cdata->evaluate(data, sel, g);
if (bDoMinMax)
static void
postprocess_item_subexpressions(t_selelem *sel)
{
+ assert(!(sel->child == NULL &&
+ (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)));
+
/* Process children. */
if (sel->type != SEL_SUBEXPRREF)
{
sel->u.cgrp.name = name;
sel->evaluate = &_gmx_sel_evaluate_subexpr_staticeval;
- if (sel->cdata)
- {
- sel->cdata->evaluate = sel->evaluate;
- }
+ sel->cdata->evaluate = sel->evaluate;
+
_gmx_selelem_free_values(sel->child);
sel->child->mempool = NULL;
_gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
char line[STRLEN];
gmx_structurefactors *gsf;
double a1,a2,a3,a4,b1,b2,b3,b4,c;
- int n,p;
+ int p;
int i;
int nralloc=10;
int line_no;
char atomn[32];
fp=libopen(datfn);
line_no = 0;
-
snew(gsf,1);
snew(gsf->atomnm,nralloc);
snew(gsf->a,nralloc);
snew(gsf->b,nralloc);
snew(gsf->c,nralloc);
- snew(gsf->n,nralloc);
snew(gsf->p,nralloc);
+ gsf->n=NULL;
gsf->nratoms=line_no;
while(get_a_line(fp,line,STRLEN)) {
i=line_no;
- if (sscanf(line,"%s %d %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
- atomn,&p,&n,&a1,&a2,&a3,&a4,&b1,&b2,&b3,&b4,&c) == 12) {
+ if (sscanf(line,"%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf",
+ atomn,&p,&a1,&a2,&a3,&a4,&b1,&b2,&b3,&b4,&c) == 11) {
gsf->atomnm[i]=strdup(atomn);
- gsf->n[i]=n;
gsf->p[i]=p;
snew(gsf->a[i],4);
snew(gsf->b[i],4);
srenew(gsf->a,nralloc);
srenew(gsf->b,nralloc);
srenew(gsf->c,nralloc);
- srenew(gsf->n,nralloc);
srenew(gsf->p,nralloc);
}
}
srenew(gsf->a,gsf->nratoms);
srenew(gsf->b,gsf->nratoms);
srenew(gsf->c,gsf->nratoms);
- srenew(gsf->n,gsf->nratoms);
srenew(gsf->p,gsf->nratoms);
fclose(fp);
return success;
}
-int atp_size (int *index_atp)
-{
- int i = 0;
-
- while (index_atp[i])
- i++;
- return i;
-}
-
-
extern int do_scattering_intensity (const char* fnTPS, const char* fnNDX,
const char* fnXVG, const char *fnTRX,
const char* fnDAT,
sfree(sf->b);
sfree(sf->atomnm);
sfree(sf->p);
- sfree(sf->n);
sfree(sf->c);
sfree(sf);
#include <ctype.h>
#include "sysstuff.h"
+#include "typedefs.h"
+#include "vmdio.h"
#include "string2.h"
#include "smalloc.h"
#include "pbc.h"
#include "confio.h"
#include "checkpoint.h"
#include "wgms.h"
-#include "vmdio.h"
#include <math.h>
/* defines for frame counter output */
for(i=0; i<nind; i++)
copy_rvec(fr->f[ind[i]],fout[i]);
}
+ /* no break */
case efXTC:
case efG87:
if (fr->bX) {
case efTRR:
if (vout) sfree(vout);
if (fout) sfree(fout);
+ /* no break */
case efXTC:
case efG87:
sfree(xout);
fr->bTime=TRUE;
fr->time=sh.t;
fr->bLambda = TRUE;
+ fr->bFepState = TRUE;
fr->lambda = sh.lambda;
fr->bBox = sh.box_size>0;
if (fr->flags & (TRX_READ_X | TRX_NEED_X)) {
return bRet;
}
-static void choose_ff(FILE *fp)
+static void choose_file_format(FILE *fp)
{
int i,m,c;
int rc;
initcount(status);
clear_mat(box);
- choose_ff(fp);
+ choose_file_format(fp);
for(m=0; (m<DIM); m++)
box[m][m]=status->BOX[m];
* accuracy of the control over -b and -e options.
*/
if (bTimeSet(TBEGIN) && (fr->time < rTimeValue(TBEGIN))) {
- if (xtc_seek_time(status->fio, rTimeValue(TBEGIN),fr->natoms)) {
- gmx_fatal(FARGS,"Specified frame doesn't exist or file not seekable");
+ if (xtc_seek_time(status->fio, rTimeValue(TBEGIN),fr->natoms,TRUE)) {
+ gmx_fatal(FARGS,"Specified frame (time %f) doesn't exist or file corrupt/inconsistent.",
+ rTimeValue(TBEGIN));
}
initcount(status);
}
bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio),fr);
break;
default:
-#ifdef GMX_DLOPEN
+#ifdef GMX_USE_PLUGINS
bRet = read_next_vmd_frame(dummy,fr);
#else
gmx_fatal(FARGS,"DEATH HORROR in read_next_frame ftp=%s,status=%s",
}
if (bRet) {
- bMissingData = ((fr->flags & TRX_NEED_X && !fr->bX) ||
- (fr->flags & TRX_NEED_V && !fr->bV) ||
- (fr->flags & TRX_NEED_F && !fr->bF));
+ bMissingData = (((fr->flags & TRX_NEED_X) && !fr->bX) ||
+ ((fr->flags & TRX_NEED_V) && !fr->bV) ||
+ ((fr->flags & TRX_NEED_F) && !fr->bF));
bSkip = FALSE;
if (!bMissingData) {
ct=check_times2(fr->time,fr->t0,fr->tpf,fr->tppf,fr->bDouble);
- if (ct == 0 || (fr->flags & TRX_DONT_SKIP && ct<0)) {
+ if (ct == 0 || ((fr->flags & TRX_DONT_SKIP) && ct<0)) {
printcount(status, oenv,fr->time,FALSE);
} else if (ct > 0)
bRet = FALSE;
bFirst = FALSE;
break;
default:
-#ifdef GMX_DLOPEN
+#ifdef GMX_USE_PLUGINS
fprintf(stderr,"The file format of %s is not a known trajectory format to GROMACS.\n"
"Please make sure that the file is a trajectory!\n"
"GROMACS will now assume it to be a trajectory and will try to open it using the VMD plug-ins.\n"
}
#else
gmx_fatal(FARGS,"Not supported in read_first_frame: %s. Please make sure that the file is a trajectory.\n"
- "GROMACS is not compiled with DLOPEN. Thus it cannot read non-GROMACS trajectory formats.\n"
- "Please compile with DLOPEN support if you want to read non-GROMACS trajectory formats.\n",fn);
+ "GROMACS is not compiled with plug-in support. Thus it cannot read non-GROMACS trajectory formats using the VMD plug-ins.\n"
+ "Please compile with plug-in support if you want to read non-GROMACS trajectory formats.\n",fn);
#endif
break;
}
clear_rvec(fr->v[i]);
}
-int read_first_v(const output_env_t oenv, t_trxstatus **status,const char *fn,
- real *t, rvec **v,matrix box)
-{
- t_trxframe fr;
-
- read_first_frame(oenv,status,fn,&fr,TRX_NEED_V);
- *t = fr.time;
- clear_v(&fr);
- *v = fr.v;
- copy_mat(fr.box,box);
-
- return fr.natoms;
-}
-
-gmx_bool read_next_v(const output_env_t oenv,t_trxstatus *status,real *t,
- int natoms,rvec v[], matrix box)
-{
- t_trxframe fr;
- gmx_bool bRet;
-
- clear_trxframe(&fr,TRUE);
- fr.flags = TRX_NEED_V;
- fr.natoms = natoms;
- fr.time = *t;
- fr.v = v;
- bRet = read_next_frame(oenv,status,&fr);
- *t = fr.time;
- clear_v(&fr);
- copy_mat(fr.box,box);
-
- return bRet;
-}
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+#include "gmx_header_config.h"
#include <stdio.h>
#include <math.h>
#include <ctype.h>
-#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
+#ifdef GMX_NATIVE_WINDOWS
#include <direct.h>
#include <io.h>
#endif
sfree(model[nwm]);
}
}
- fclose(fp);
+ ffclose(fp);
fprintf(stderr,"%2d: %s\n",nwm+1,"None");
do
else
{
strncpy(ffdir_parent,ffdir,STRLEN-1);
+ ffdir_parent[STRLEN-1]='\0'; /*make sure it is 0-terminated even for long string*/
p=strrchr(ffdir_parent,'/');
*p='\0';
}
}
-static gmx_bool inter_res_bond(const t_rbonded *b)
-{
- return (b->AI[0] == '-' || b->AI[0] == '+' ||
- b->AJ[0] == '-' || b->AJ[0] == '+');
-}
-
static void at2bonds(t_params *psb, t_hackblock *hb,
t_atoms *atoms,
rvec x[],
if (bRM && ((tern >= 0 && ntdb[tern] == NULL) ||
(terc >= 0 && ctdb[terc] == NULL))) {
- gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Edit a .n.tdb and/or .c.tdb file.");
+ gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends and the force field does not provide terminal entries or files. Fix your terminal residues so that they match the residue database (.rtp) entries, or provide terminal database entries (.tdb).");
}
if (bRM && ((tern >= 0 && ntdb[tern]->nhack == 0) ||
(terc >= 0 && ctdb[terc]->nhack == 0))) {
- gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends. Select a proper terminal entry.");
+ gmx_fatal(FARGS,"There is a dangling bond at at least one of the terminal ends. Fix your coordinate file, add a new terminal database entry (.tdb), or select the proper existing terminal entry.");
}
}
int nres = atoms->nres;
gmx_bool bAddCMAP;
atom_id cmap_atomid[NUM_CMAP_ATOMS];
- int cmap_chainnum, this_residue_index;
+ int cmap_chainnum=-1, this_residue_index;
if (debug)
ptr = "cmap";
#endif
#include <stdio.h>
+#include <assert.h>
#include "domdec.h"
#include "network.h"
#include "perf_est.h"
return -1;
}
+ assert(ddbox->npbcdim<=DIM);
+
/* Check if the triclinic requirements are met */
for(i=0; i<DIM; i++)
{
{
if (cr->npmenodes < 0)
{
- if (cr->nnodes <= 10)
+ /* Use PME nodes when the number of nodes is more than 16 */
+ if (cr->nnodes <= 18)
{
cr->npmenodes = 0;
+ if (fplog)
+ {
+ fprintf(fplog,"Using %d separate PME nodes, as there are too few total\n nodes for efficient splitting\n",cr->npmenodes);
+ }
}
else
{
cr->npmenodes = guess_npme(fplog,mtop,ir,box,cr->nnodes);
+ if (fplog)
+ {
+ fprintf(fplog,"Using %d separate PME nodes, as guessed by mdrun\n",cr->npmenodes);
+ }
}
}
- if (fplog)
+ else
{
- fprintf(fplog,"Using %d separate PME nodes\n",cr->npmenodes);
+ if (fplog)
+ {
+ fprintf(fplog,"Using %d separate PME nodes, per user request\n",cr->npmenodes);
+ }
}
}
/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
- *
+ *
* If you want to redistribute modifications, please consider that
* scientific software is very special. Version control is crucial -
* bugs must be traceable. We will be happy to consider code for
* inclusion in the official distribution, but derived work must not
* be called official GROMACS. Details are found in the README & COPYING
* files - if they are missing, get the official version at www.gromacs.org.
- *
+ *
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
- *
+ *
* For more info, check our website at http://www.gromacs.org
- *
+ *
* And Hey:
* GROwing Monsters And Cloning Shrimps
*/
#include "mtop_util.h"
#include "xvgr.h"
#include "gmxfio.h"
-
+#include "mdrun.h"
#include "mdebin_bar.h"
static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
-static const char *tricl_boxs_nm[] = {
+static const char *tricl_boxs_nm[] = {
"Box-XX", "Box-YY", "Box-ZZ",
- "Box-YX", "Box-ZX", "Box-ZY"
+ "Box-YX", "Box-ZX", "Box-ZY"
};
static const char *vol_nm[] = { "Volume" };
#define NBOXS asize(boxs_nm)
#define NTRICLBOXS asize(tricl_boxs_nm)
-
t_mdebin *init_mdebin(ener_file_t fp_ene,
const gmx_mtop_t *mtop,
const t_inputrec *ir,
snew(md,1);
+ md->bVir=TRUE;
+ md->bPress=TRUE;
+ md->bSurft=TRUE;
+ md->bMu=TRUE;
+
if (EI_DYNAMICS(ir->eI))
{
md->delta_t = ir->delta_t;
md->bEner[i] = b14;
else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
md->bEner[i] = FALSE;
- else if ((i == F_DVDL) || (i == F_DKDL))
+ else if ((i == F_DVDL_COUL && ir->fepvals->separate_dvdl[efptCOUL]) ||
+ (i == F_DVDL_VDW && ir->fepvals->separate_dvdl[efptVDW]) ||
+ (i == F_DVDL_BONDED && ir->fepvals->separate_dvdl[efptBONDED]) ||
+ (i == F_DVDL_RESTRAINT && ir->fepvals->separate_dvdl[efptRESTRAINT]) ||
+ (i == F_DKDL && ir->fepvals->separate_dvdl[efptMASS]) ||
+ (i == F_DVDL && ir->fepvals->separate_dvdl[efptFEP]))
md->bEner[i] = (ir->efep != efepNO);
- else if (i == F_DHDL_CON)
- md->bEner[i] = (ir->efep != efepNO && md->bConstr);
else if ((interaction_function[i].flags & IF_VSITE) ||
(i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
md->bEner[i] = FALSE;
md->bEner[i] = FALSE;
else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
md->bEner[i] = EI_DYNAMICS(ir->eI);
- else if (i==F_VTEMP)
+ else if (i==F_VTEMP)
md->bEner[i] = (EI_DYNAMICS(ir->eI) && getenv("GMX_VIRIAL_TEMPERATURE"));
else if (i == F_DISPCORR || i == F_PDISPCORR)
md->bEner[i] = (ir->eDispCorr != edispcNO);
else if (i == F_CONNBONDS)
md->bEner[i] = FALSE;
else if (i == F_COM_PULL)
- md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F);
+ md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot);
else if (i == F_ECONSERVED)
md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
(ir->epc == epcNO || ir->epc==epcMTTK));
md->bEner[F_TEMP] = TRUE;
#endif
+ /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
+ if (ir->bAdress && !debug) {
+ for (i = 0; i < F_NRE; i++) {
+ md->bEner[i] = FALSE;
+ if(i == F_EKIN){ md->bEner[i] = TRUE;}
+ if(i == F_TEMP){ md->bEner[i] = TRUE;}
+ }
+ md->bVir=FALSE;
+ md->bPress=FALSE;
+ md->bSurft=FALSE;
+ md->bMu=FALSE;
+ }
+
md->f_nre=0;
for(i=0; i<F_NRE; i++)
{
md->bDynBox = DYNAMIC_BOX(*ir);
md->etc = ir->etc;
md->bNHC_trotter = IR_NVT_TROTTER(ir);
- md->bMTTK = IR_NPT_TROTTER(ir);
+ md->bPrintNHChains = ir-> bPrintNHChains;
+ md->bMTTK = (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir));
md->ebin = mk_ebin();
/* Pass NULL for unit to let get_ebin_space determine the units
}
if (md->bDynBox)
{
- md->ib = get_ebin_space(md->ebin,
- md->bTricl ? NTRICLBOXS : NBOXS,
+ md->ib = get_ebin_space(md->ebin,
+ md->bTricl ? NTRICLBOXS : NBOXS,
md->bTricl ? tricl_boxs_nm : boxs_nm,
unit_length);
md->ivol = get_ebin_space(md->ebin, 1, vol_nm, unit_volume);
md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
}
- md->ivir = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
- md->ipres = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
- md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
+ if (md->bVir)
+ md->ivir = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
+ if (md->bPress)
+ md->ipres = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
+ if (md->bSurft)
+ md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
unit_surft_bar);
if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
{
md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
boxvel_nm,unit_vel);
}
- md->imu = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
+ if (md->bMu)
+ md->imu = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
if (ir->cos_accel != 0)
{
md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
}
n=groups->grps[egcENER].nr;
- md->nEg=n;
- md->nE=(n*(n+1))/2;
+ /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
+ if (!ir->bAdress){
+ /*standard simulation*/
+ md->nEg=n;
+ md->nE=(n*(n+1))/2;
+ }
+ else if (!debug) {
+ /*AdResS simulation*/
+ md->nU=0;
+ md->nEg=0;
+ md->nE=0;
+ md->nEc=0;
+ md->isvir=FALSE;
+ }
snew(md->igrp,md->nE);
if (md->nE > 1)
{
}
md->nTC=groups->grps[egcTC].nr;
- md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
+ md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
if (md->bMTTK)
{
- md->nTCP = 1; /* assume only one possible coupling system for barostat
+ md->nTCP = 1; /* assume only one possible coupling system for barostat
for now */
- }
- else
+ }
+ else
{
md->nTCP = 0;
}
-
- if (md->etc == etcNOSEHOOVER) {
- if (md->bNHC_trotter) {
+ if (md->etc == etcNOSEHOOVER)
+ {
+ if (md->bNHC_trotter)
+ {
md->mde_n = 2*md->nNHC*md->nTC;
}
- else
+ else
{
md->mde_n = 2*md->nTC;
}
{
md->mdeb_n = 2*md->nNHC*md->nTCP;
}
- } else {
+ } else {
md->mde_n = md->nTC;
md->mdeb_n = 0;
}
md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
unit_temp_K);
- bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
-
if (md->etc == etcNOSEHOOVER)
{
- if (bNoseHoover)
+ if (md->bPrintNHChains)
{
- if (md->bNHC_trotter)
+ if (md->bNHC_trotter)
{
- for(i=0; (i<md->nTC); i++)
+ for(i=0; (i<md->nTC); i++)
{
ni=groups->grps[egcTC].nm_ind[i];
bufi = *(groups->grpname[ni]);
- for(j=0; (j<md->nNHC); j++)
+ for(j=0; (j<md->nNHC); j++)
{
sprintf(buf,"Xi-%d-%s",j,bufi);
grpnms[2*(i*md->nNHC+j)]=strdup(buf);
}
md->itc=get_ebin_space(md->ebin,md->mde_n,
(const char **)grpnms,unit_invtime);
- if (md->bMTTK)
+ if (md->bMTTK)
{
- for(i=0; (i<md->nTCP); i++)
+ for(i=0; (i<md->nTCP); i++)
{
bufi = baro_nm[0]; /* All barostat DOF's together for now. */
- for(j=0; (j<md->nNHC); j++)
+ for(j=0; (j<md->nNHC); j++)
{
sprintf(buf,"Xi-%d-%s",j,bufi);
grpnms[2*(i*md->nNHC+j)]=strdup(buf);
md->itcb=get_ebin_space(md->ebin,md->mdeb_n,
(const char **)grpnms,unit_invtime);
}
- }
+ }
else
{
- for(i=0; (i<md->nTC); i++)
+ for(i=0; (i<md->nTC); i++)
{
ni=groups->grps[egcTC].nm_ind[i];
bufi = *(groups->grpname[ni]);
}
}
}
- else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
+ else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
md->etc == etcVRESCALE)
{
for(i=0; (i<md->nTC); i++)
md->print_grpnms=NULL;
/* check whether we're going to write dh histograms */
- md->dhc=NULL;
- if (ir->separate_dhdl_file == sepdhdlfileNO )
+ md->dhc=NULL;
+ if (ir->fepvals->separate_dhdl_file == esepdhdlfileNO )
{
- snew(md->dhc, 1);
+ /* Currently dh histograms are only written with dynamics */
+ if (EI_DYNAMICS(ir->eI))
+ {
+ snew(md->dhc, 1);
- mde_delta_h_coll_init(md->dhc, ir);
+ mde_delta_h_coll_init(md->dhc, ir);
+ }
md->fp_dhdl = NULL;
}
else
{
md->fp_dhdl = fp_dhdl;
}
- md->dhdl_derivatives = (ir->dhdl_derivatives==dhdlderivativesYES);
+ if (ir->bSimTemp) {
+ int i;
+ snew(md->temperatures,ir->fepvals->n_lambda);
+ for (i=0;i<ir->fepvals->n_lambda;i++)
+ {
+ md->temperatures[i] = ir->simtempvals->temperatures[i];
+ }
+ }
return md;
}
-FILE *open_dhdl(const char *filename,const t_inputrec *ir,
- const output_env_t oenv)
+extern FILE *open_dhdl(const char *filename,const t_inputrec *ir,
+ const output_env_t oenv)
{
FILE *fp;
- const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
+ const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda",
+ *lambdastate="\\lambda state",*remain="remaining";
char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
+ int i,np,nps,nsets,nsets_de,nsetsbegin;
+ t_lambda *fep;
char **setname;
char buf[STRLEN];
+ int bufplace=0;
+
+ int nsets_dhdl = 0;
+ int s = 0;
+ int nsetsextend;
- sprintf(label_x,"%s (%s)","Time",unit_time);
- if (ir->n_flambda == 0)
+ /* for simplicity */
+ fep = ir->fepvals;
+
+ if (fep->n_lambda == 0)
{
sprintf(title,"%s",dhdl);
+ sprintf(label_x,"Time (ps)");
sprintf(label_y,"%s (%s %s)",
dhdl,unit_energy,"[\\lambda]\\S-1\\N");
}
else
{
- sprintf(title,"%s, %s",dhdl,deltag);
- sprintf(label_y,"(%s)",unit_energy);
+ sprintf(title,"%s and %s",dhdl,deltag);
+ sprintf(label_x,"Time (ps)");
+ sprintf(label_y,"%s and %s (%s %s)",
+ dhdl,deltag,unit_energy,"[\\8l\\4]\\S-1\\N");
}
fp = gmx_fio_fopen(filename,"w+");
xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
- if (ir->delta_lambda == 0)
+ if (!(ir->bSimTemp))
{
- sprintf(buf,"T = %g (K), %s = %g",
- ir->opts.ref_t[0],lambda,ir->init_lambda);
+ bufplace = sprintf(buf,"T = %g (K) ",
+ ir->opts.ref_t[0]);
}
- else
+ if (ir->efep != efepSLOWGROWTH)
{
- sprintf(buf,"T = %g (K)",
- ir->opts.ref_t[0]);
+ if (fep->n_lambda == 0)
+ {
+ sprintf(&(buf[bufplace]),"%s = %g",
+ lambda,fep->init_lambda);
+ }
+ else
+ {
+ sprintf(&(buf[bufplace]),"%s = %d",
+ lambdastate,fep->init_fep_state);
+ }
}
xvgr_subtitle(fp,buf,oenv);
- if (ir->n_flambda > 0)
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i]) {nsets_dhdl++;}
+ }
+
+ /* count the number of delta_g states */
+ nsets_de = fep->n_lambda;
+
+ nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */
+
+ if (fep->n_lambda>0 && ir->bExpanded)
+ {
+ nsets += 1; /*add fep state for expanded ensemble */
+ }
+
+ if (fep->bPrintEnergy)
+ {
+ nsets += 1; /* add energy to the dhdl as well */
+ }
+
+ nsetsextend = nsets;
+ if ((ir->epc!=epcNO) && (fep->n_lambda>0))
+ {
+ nsetsextend += 1; /* for PV term, other terms possible if required for the reduced potential (only needed with foreign lambda) */
+ }
+ snew(setname,nsetsextend);
+
+ if (ir->bExpanded)
+ {
+ /* state for the fep_vals, if we have alchemical sampling */
+ sprintf(buf,"%s","Thermodynamic state");
+ setname[s] = strdup(buf);
+ s+=1;
+ }
+
+ if (fep->bPrintEnergy)
+ {
+ sprintf(buf,"%s (%s)","Energy",unit_energy);
+ setname[s] = strdup(buf);
+ s+=1;
+ }
+
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i]) {
+ sprintf(buf,"%s (%s)",dhdl,efpt_names[i]);
+ setname[s] = strdup(buf);
+ s+=1;
+ }
+ }
+
+ if (fep->n_lambda > 0)
{
- int nsets,s,nsi=0;
/* g_bar has to determine the lambda values used in this simulation
- * from this xvg legend. */
- nsets = ( (ir->dhdl_derivatives==dhdlderivativesYES) ? 1 : 0) +
- ir->n_flambda;
- snew(setname,nsets);
- if (ir->dhdl_derivatives == dhdlderivativesYES)
+ * from this xvg legend.
+ */
+
+ if (ir->bExpanded) {
+ nsetsbegin = 1; /* for including the expanded ensemble */
+ } else {
+ nsetsbegin = 0;
+ }
+
+ if (fep->bPrintEnergy)
{
- sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
- setname[nsi++] = gmx_strdup(buf);
+ nsetsbegin += 1;
}
- for(s=0; s<ir->n_flambda; s++)
+ nsetsbegin += nsets_dhdl;
+
+ for(s=nsetsbegin; s<nsets; s++)
{
- sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s]);
- setname[nsi++] = gmx_strdup(buf);
+ nps = sprintf(buf,"%s %s (",deltag,lambda);
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ np = sprintf(&buf[nps],"%g,",fep->all_lambda[i][s-(nsetsbegin)]);
+ nps += np;
+ }
+ }
+ if (ir->bSimTemp)
+ {
+ /* print the temperature for this state if doing simulated annealing */
+ sprintf(&buf[nps],"T = %g (%s))",ir->simtempvals->temperatures[s-(nsetsbegin)],unit_temp_K);
+ }
+ else
+ {
+ sprintf(&buf[nps-1],")"); /* -1 to overwrite the last comma */
+ }
+ setname[s] = strdup(buf);
+ }
+ if (ir->epc!=epcNO) {
+ np = sprintf(buf,"pV (%s)",unit_energy);
+ setname[nsetsextend-1] = strdup(buf); /* the first entry after nsets */
}
- xvgr_legend(fp,nsets,(const char**)setname,oenv);
- for(s=0; s<nsets; s++)
+ xvgr_legend(fp,nsetsextend,(const char **)setname,oenv);
+
+ for(s=0; s<nsetsextend; s++)
{
sfree(setname[s]);
}
for(i=j=0; (i<F_NRE); i++)
if (md->bEner[i])
ecpy[j++] = e[i];
- if (j != md->f_nre)
+ if (j != md->f_nre)
gmx_incons("Number of energy terms wrong");
}
-void upd_mdebin(t_mdebin *md, gmx_bool write_dhdl,
+void upd_mdebin(t_mdebin *md,
+ gmx_bool bDoDHDL,
gmx_bool bSum,
double time,
real tmass,
gmx_enerdata_t *enerd,
t_state *state,
+ t_lambda *fep,
+ t_expanded *expand,
matrix box,
tensor svir,
tensor fvir,
real bs[NTRICLBOXS],vol,dens,pv,enthalpy;
real eee[egNR];
real ecopy[F_NRE];
+ double store_dhdl[efptNR];
+ double *dE=NULL;
+ real store_energy=0;
real tmp;
- gmx_bool bNoseHoover;
/* Do NOT use the box in the state variable, but the separate box provided
* as an argument. This is because we sometimes need to write the box from
}
vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
-
add_ebin(md->ebin,md->ib ,nboxs,bs ,bSum);
add_ebin(md->ebin,md->ivol ,1 ,&vol ,bSum);
add_ebin(md->ebin,md->idens,1 ,&dens,bSum);
if (md->bDiagPres)
{
/* This is pV (in kJ/mol). The pressure is the reference pressure,
- not the instantaneous pressure */
+ not the instantaneous pressure */
pv = vol*md->ref_p/PRESFAC;
add_ebin(md->ebin,md->ipv ,1 ,&pv ,bSum);
add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
}
- add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
- add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
- tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
- add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
+ if (md->bVir)
+ add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
+ if (md->bPress)
+ add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
+ if (md->bSurft){
+ tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
+ add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
+ }
if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
{
tmp6[0] = state->boxv[XX][XX];
tmp6[5] = state->boxv[ZZ][YY];
add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
}
- add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
+ if(md->bMu)
+ add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
if (ekind && ekind->cosacc.cos_accel != 0)
{
vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
/* 1/viscosity, unit 1/(kg m^-1 s^-1) */
tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
*dens*vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
- add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
+ add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
}
if (md->nE > 1)
{
}
add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
- /* whether to print Nose-Hoover chains: */
- bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL);
-
if (md->etc == etcNOSEHOOVER)
{
- if (bNoseHoover)
+ /* whether to print Nose-Hoover chains: */
+ if (md->bPrintNHChains)
{
if (md->bNHC_trotter)
{
- for(i=0; (i<md->nTC); i++)
+ for(i=0; (i<md->nTC); i++)
{
- for (j=0;j<md->nNHC;j++)
+ for (j=0;j<md->nNHC;j++)
{
k = i*md->nNHC+j;
md->tmp_r[2*k] = state->nosehoover_xi[k];
md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
}
}
- add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
+ add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
if (md->bMTTK) {
- for(i=0; (i<md->nTCP); i++)
+ for(i=0; (i<md->nTCP); i++)
{
- for (j=0;j<md->nNHC;j++)
+ for (j=0;j<md->nNHC;j++)
{
k = i*md->nNHC+j;
md->tmp_r[2*k] = state->nhpres_xi[k];
md->tmp_r[2*k+1] = state->nhpres_vxi[k];
}
}
- add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
+ add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
}
-
- }
- else
+ }
+ else
{
for(i=0; (i<md->nTC); i++)
{
}
}
}
- else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
+ else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
md->etc == etcVRESCALE)
{
for(i=0; (i<md->nTC); i++)
ebin_increase_count(md->ebin,bSum);
/* BAR + thermodynamic integration values */
- if (write_dhdl)
+ if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda > 0))
{
- if (md->fp_dhdl)
- {
- fprintf(md->fp_dhdl,"%.4f", time);
-
- if (md->dhdl_derivatives)
+ snew(dE,enerd->n_lambda-1);
+ for(i=0; i<enerd->n_lambda-1; i++) {
+ dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0]; /* zero for simulated tempering */
+ if (md->temperatures!=NULL)
{
- fprintf(md->fp_dhdl," %g", enerd->term[F_DVDL]+
- enerd->term[F_DKDL]+
- enerd->term[F_DHDL_CON]);
+ /* MRS: is this right, given the way we have defined the exchange probabilities? */
+ /* is this even useful to have at all? */
+ dE[i] += (md->temperatures[i]/md->temperatures[state->fep_state]-1.0)*enerd->term[F_EKIN];
}
- for(i=1; i<enerd->n_lambda; i++)
+ }
+ }
+
+ if (md->fp_dhdl && bDoDHDL)
+ {
+ fprintf(md->fp_dhdl,"%.4f",time);
+ /* the current free energy state */
+
+ /* print the current state if we are doing expanded ensemble */
+ if (expand->elmcmove > elmcmoveNO) {
+ fprintf(md->fp_dhdl," %4d",state->fep_state);
+ }
+ /* total energy (for if the temperature changes */
+ if (fep->bPrintEnergy)
+ {
+ store_energy = enerd->term[F_ETOT];
+ fprintf(md->fp_dhdl," %#.8g",store_energy);
+ }
+
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
{
- fprintf(md->fp_dhdl," %g",
- enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
+ fprintf(md->fp_dhdl," %#.8g",enerd->term[F_DVDL+i]); /* assumes F_DVDL is first */
}
- fprintf(md->fp_dhdl,"\n");
}
- /* and the binary BAR output */
- if (md->dhc)
+ for(i=1; i<enerd->n_lambda; i++)
{
- mde_delta_h_coll_add_dh(md->dhc,
- enerd->term[F_DVDL]+ enerd->term[F_DKDL]+
- enerd->term[F_DHDL_CON],
- enerd->enerpart_lambda, time,
- state->lambda);
+ fprintf(md->fp_dhdl," %#.8g",dE[i-1]);
+
+ }
+ if ((md->epc!=epcNO) && (enerd->n_lambda > 0))
+ {
+ fprintf(md->fp_dhdl," %#.8g",pv); /* PV term only needed when there are alternate state lambda */
+ }
+ fprintf(md->fp_dhdl,"\n");
+ /* and the binary free energy output */
+ }
+ if (md->dhc && bDoDHDL)
+ {
+ int idhdl = 0;
+ for (i=0;i<efptNR;i++)
+ {
+ if (fep->separate_dvdl[i])
+ {
+ store_dhdl[idhdl] = enerd->term[F_DVDL+i]; /* assumes F_DVDL is first */
+ idhdl+=1;
+ }
}
+ /* store_dh is dE */
+ mde_delta_h_coll_add_dh(md->dhc,
+ (double)state->fep_state,
+ store_energy,
+ pv,
+ (expand->elamstats>elamstatsNO),
+ (fep->bPrintEnergy),
+ (md->epc!=epcNO),
+ idhdl,
+ fep->n_lambda,
+ store_dhdl,
+ dE,
+ time);
+ }
+ if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda >0))
+ {
+ sfree(dE);
}
}
+
void upd_mdebin_step(t_mdebin *md)
{
- ebin_increase_count(md->ebin,FALSE);
+ ebin_increase_count(md->ebin,FALSE);
}
static void npr(FILE *log,int n,char c)
fprintf(log,"\n");
}
-void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
+void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lambda)
{
char buf[22];
fprintf(log," %12s %12s %12s\n"
" %12s %12.5f %12.5f\n\n",
- "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
+ "Step","Time","Lambda",gmx_step_str(steps,buf),time,lambda);
}
void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
real *block[enxNR];
/* temporary arrays for the lambda values to write out */
- double enxlambda_data[2];
+ double enxlambda_data[2];
t_enxframe fr;
nr[enxOR] = fcd->orires.nr;
block[enxOR] = fcd->orires.otav;
id[enxOR] = enxOR;
- nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
+ nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
fcd->orires.nr : 0;
block[enxORI] = fcd->orires.oinsl;
id[enxORI] = enxORI;
nr[enxORT] = fcd->orires.nex*12;
block[enxORT] = fcd->orires.eig;
id[enxORT] = enxORT;
- }
+ }
/* whether we are going to wrte anything out: */
if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
for(b=0;b<fr.nblock;b++)
{
add_subblocks_enxblock(&(fr.block[b]), 1);
- fr.block[b].id=id[b];
+ fr.block[b].id=id[b];
fr.block[b].sub[0].nr = nr[b];
#ifndef GMX_DOUBLE
fr.block[b].sub[0].type = xdr_datatype_float;
mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
}
+ /* we can now free & reset the data in the blocks */
+ if (md->dhc)
+ {
+ mde_delta_h_coll_reset(md->dhc);
+ }
+
/* do the actual I/O */
do_enx(fp_ene,&fr);
gmx_fio_check_file_position(enx_file_pointer(fp_ene));
/* We have stored the sums, so reset the sum history */
reset_ebin_sums(md->ebin);
}
-
- /* we can now free & reset the data in the blocks */
- if (md->dhc)
- mde_delta_h_coll_reset(md->dhc);
}
free_enxframe(&fr);
break;
print_orires_log(log,&(fcd->orires));
}
fprintf(log," Energies (%s)\n",unit_energy);
- pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
+ pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
fprintf(log,"\n");
if (!bCompact)
if (md->bDynBox)
{
pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
- mode,TRUE);
+ mode,TRUE);
fprintf(log,"\n");
}
if (md->bConstrVir)
{
fprintf(log," Constraint Virial (%s)\n",unit_energy);
- pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
+ pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
fprintf(log,"\n");
fprintf(log," Force Virial (%s)\n",unit_energy);
- pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
+ pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ }
+ if (md->bVir)
+ {
+ fprintf(log," Total Virial (%s)\n",unit_energy);
+ pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
+ fprintf(log,"\n");
+ }
+ if (md->bPress)
+ {
+ fprintf(log," Pressure (%s)\n",unit_pres_bar);
+ pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
fprintf(log,"\n");
}
- fprintf(log," Total Virial (%s)\n",unit_energy);
- pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
- fprintf(log,"\n");
- fprintf(log," Pressure (%s)\n",unit_pres_bar);
- pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
- fprintf(log,"\n");
fprintf(log," Total Dipole (%s)\n",unit_dipole_D);
- pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
+ pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
fprintf(log,"\n");
if (md->nE > 1)
(enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
}
if (mdebin->dhc)
- {
+ {
mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
}
}
/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
*
- *
+ *
* This source code is part of
- *
+ *
* G R O M A C S
- *
+ *
* GROningen MAchine for Chemical Simulations
- *
+ *
* VERSION 3.2.0
* Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
- *
+ *
* If you want to redistribute modifications, please consider that
* scientific software is very special. Version control is crucial -
* bugs must be traceable. We will be happy to consider code for
* inclusion in the official distribution, but derived work must not
* be called official GROMACS. Details are found in the README & COPYING
* files - if they are missing, get the official version at www.gromacs.org.
- *
+ *
* To help us fund GROMACS development, we humbly ask that you cite
* the papers on the package - you can find them in the top README file.
- *
+ *
* For more info, check our website at http://www.gromacs.org
- *
+ *
* And Hey:
* GROwing Monsters And Cloning Shrimps
*/
static em_state_t *init_em_state()
{
em_state_t *ems;
-
+
snew(ems,1);
+ /* does this need to be here? Should the array be declared differently (staticaly)in the state definition? */
+ snew(ems->s.lambda,efptNR);
+
return ems;
}
if (bDone)
fprintf(fp,"\n%s converged to Fmax < %g in %s steps\n",
- alg,ftol,gmx_step_str(count,buf));
+ alg,ftol,gmx_step_str(count,buf));
else if(count<nsteps)
fprintf(fp,"\n%s converged to machine precision in %s steps,\n"
"but did not reach the requested Fmax < %g.\n",
alg,gmx_step_str(count,buf),ftol);
- else
+ else
fprintf(fp,"\n%s did not converge to Fmax < %g in %s steps.\n",
alg,ftol,gmx_step_str(count,buf));
#ifdef GMX_DOUBLE
- fprintf(fp,"Potential Energy = %21.14e\n",epot);
- fprintf(fp,"Maximum force = %21.14e on atom %d\n",fmax,nfmax+1);
- fprintf(fp,"Norm of force = %21.14e\n",fnorm);
+ fprintf(fp,"Potential Energy = %21.14e\n",epot);
+ fprintf(fp,"Maximum force = %21.14e on atom %d\n",fmax,nfmax+1);
+ fprintf(fp,"Norm of force = %21.14e\n",fnorm);
#else
- fprintf(fp,"Potential Energy = %14.7e\n",epot);
- fprintf(fp,"Maximum force = %14.7e on atom %d\n",fmax,nfmax+1);
- fprintf(fp,"Norm of force = %14.7e\n",fnorm);
+ fprintf(fp,"Potential Energy = %14.7e\n",epot);
+ fprintf(fp,"Maximum force = %14.7e on atom %d\n",fmax,nfmax+1);
+ fprintf(fp,"Norm of force = %14.7e\n",fnorm);
#endif
}
{
int start,homenr,i;
real dvdlambda;
-
+
if (fplog)
{
fprintf(fplog,"Initiating %s\n",title);
}
-
+
state_global->ngtc = 0;
-
- /* Initiate some variables */
- if (ir->efep != efepNO)
- {
- state_global->lambda = ir->init_lambda;
- }
- else
- {
- state_global->lambda = 0.0;
- }
-
+
+ /* Initialize lambda variables */
+ initialize_lambdas(fplog,ir,&(state_global->fep_state),state_global->lambda,NULL);
+
init_nrnb(nrnb);
-
+
if (DOMAINDECOMP(cr))
{
*top = dd_init_local_top(top_global);
-
+
dd_init_local_state(cr->dd,state_global,&ems->s);
*f = NULL;
-
+
/* Distribute the charge groups over the nodes from the master node */
dd_partition_system(fplog,ir->init_step,cr,TRUE,1,
state_global,top_global,ir,
fr,vsite,NULL,constr,
nrnb,NULL,FALSE);
dd_store_state(cr->dd,&ems->s);
-
+
if (ir->nstfout)
{
snew(*f_global,top_global->natoms);
copy_rvec(state_global->x[i],ems->s.x[i]);
}
copy_mat(state_global->box,ems->s.box);
-
+
if (PAR(cr) && ir->eI != eiNM)
{
/* Initialize the particle decomposition and split the topology */
*top = split_system(fplog,top_global,ir,cr);
-
+
pd_cg_range(cr,&fr->cg0,&fr->hcg);
}
else
*top = gmx_mtop_generate_local_top(top_global,ir);
}
*f_global = *f;
-
+
if (ir->ePBC != epbcNONE && !ir->bPeriodicMols)
{
*graph = mk_graph(fplog,&((*top)->idef),0,top_global->natoms,FALSE,FALSE);
homenr = top_global->natoms;
}
atoms2md(top_global,ir,0,NULL,start,homenr,mdatoms);
- update_mdatoms(mdatoms,state_global->lambda);
-
+ update_mdatoms(mdatoms,state_global->lambda[efptFEP]);
+
if (vsite)
{
set_vsite_top(vsite,*top,mdatoms,cr);
}
}
-
+
if (constr)
{
if (ir->eConstrAlg == econtSHAKE &&
gmx_fatal(FARGS,"Can not do energy minimization with %s, use %s\n",
econstr_names[econtSHAKE],econstr_names[econtLINCS]);
}
-
+
if (!DOMAINDECOMP(cr))
{
set_constraints(constr,*top,ir,mdatoms,cr);
constrain(PAR(cr) ? NULL : fplog,TRUE,TRUE,constr,&(*top)->idef,
ir,NULL,cr,-1,0,mdatoms,
ems->s.x,ems->s.x,NULL,ems->s.box,
- ems->s.lambda,&dvdlambda,
+ ems->s.lambda[efptFEP],&dvdlambda,
NULL,NULL,nrnb,econqCoord,FALSE,0,0);
}
}
-
+
if (PAR(cr))
{
*gstat = global_stat_init(ir);
}
-
+
*outf = init_mdoutf(nfile,fnm,0,cr,ir,NULL);
snew(*enerd,1);
- init_enerdata(top_global->groups.grps[egcENER].nr,ir->n_flambda,*enerd);
+ init_enerdata(top_global->groups.grps[egcENER].nr,ir->fepvals->n_lambda,
+ *enerd);
if (mdebin != NULL)
{
/* Init bin for energy stuff */
- *mdebin = init_mdebin((*outf)->fp_ene,top_global,ir,NULL);
+ *mdebin = init_mdebin((*outf)->fp_ene,top_global,ir,NULL);
}
clear_rvec(mu_tot);
copy_em_coords(state,state_global);
f_global = state->f;
}
-
+
mdof_flags = 0;
if (bX) { mdof_flags |= MDOF_X; }
if (bF) { mdof_flags |= MDOF_F; }
write_traj(fplog,cr,outf,mdof_flags,
top_global,step,(double)step,
&state->s,state_global,state->f,f_global,NULL,NULL);
-
+
if (confout != NULL && MASTER(cr))
{
if (ir->ePBC != epbcNONE && !ir->bPeriodicMols && DOMAINDECOMP(cr))
if (s2->flags & (1<<estCGP))
srenew(s2->cg_p, s1->nalloc);
}
-
+
s2->natoms = s1->natoms;
- s2->lambda = s1->lambda;
+ /* Copy free energy state -> is this necessary? */
+ for (i=0;i<efptNR;i++)
+ {
+ s2->lambda[i] = s1->lambda[i];
+ }
copy_mat(s1->box,s2->box);
start = md->start;
if (constr) {
wallcycle_start(wcycle,ewcCONSTR);
dvdlambda = 0;
- constrain(NULL,TRUE,TRUE,constr,&top->idef,
+ constrain(NULL,TRUE,TRUE,constr,&top->idef,
ir,NULL,cr,count,0,md,
- s1->x,s2->x,NULL,s2->box,s2->lambda,
+ s1->x,s2->x,NULL,s2->box,s2->lambda[efptBONDED],
&dvdlambda,NULL,NULL,nrnb,econqCoord,FALSE,0,0);
wallcycle_stop(wcycle,ewcCONSTR);
}
}
-static void do_x_step(t_commrec *cr,int n,rvec *x1,real a,rvec *f,rvec *x2)
-
-{
- int start,end,i,m;
-
- if (DOMAINDECOMP(cr)) {
- start = 0;
- end = cr->dd->nat_home;
- } else if (PARTDECOMP(cr)) {
- pd_at_range(cr,&start,&end);
- } else {
- start = 0;
- end = n;
- }
-
- for(i=start; i<end; i++) {
- for(m=0; m<DIM; m++) {
- x2[i][m] = x1[i][m] + a*f[i][m];
- }
- }
-}
-
-static void do_x_sub(t_commrec *cr,int n,rvec *x1,rvec *x2,real a,rvec *f)
-
-{
- int start,end,i,m;
-
- if (DOMAINDECOMP(cr)) {
- start = 0;
- end = cr->dd->nat_home;
- } else if (PARTDECOMP(cr)) {
- pd_at_range(cr,&start,&end);
- } else {
- start = 0;
- end = n;
- }
-
- for(i=start; i<end; i++) {
- for(m=0; m<DIM; m++) {
- f[i][m] = (x1[i][m] - x2[i][m])*a;
- }
- }
-}
-
static void em_dd_partition_system(FILE *fplog,int step,t_commrec *cr,
gmx_mtop_t *top_global,t_inputrec *ir,
em_state_t *ems,gmx_localtop_t *top,
dd_store_state(cr->dd,&ems->s);
wallcycle_stop(wcycle,ewcDOMDEC);
}
-
+
static void evaluate_energy(FILE *fplog,gmx_bool bVerbose,t_commrec *cr,
t_state *state_global,gmx_mtop_t *top_global,
em_state_t *ems,gmx_localtop_t *top,
gmx_bool bNS;
int nabnsb;
tensor force_vir,shake_vir,ekin;
- real dvdl,prescorr,enercorr,dvdlcorr;
+ real dvdlambda,prescorr,enercorr,dvdlcorr;
real terminate=0;
-
+
/* Set the time to the initial time, the time does not change during EM */
t = inputrec->init_t;
nrnb,wcycle);
}
}
-
+
/* Calc force & energy on new trial position */
/* do_force always puts the charge groups in the box and shifts again
* We do not unshift, so molecules are always whole in congrad.c
ems->s.lambda,graph,fr,vsite,mu_tot,t,NULL,NULL,TRUE,
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES | GMX_FORCE_VIRIAL |
(bNS ? GMX_FORCE_NS | GMX_FORCE_DOLR : 0));
-
- /* Clear the unused shake virial and pressure */
- clear_mat(shake_vir);
- clear_mat(pres);
+
+ /* Clear the unused shake virial and pressure */
+ clear_mat(shake_vir);
+ clear_mat(pres);
/* Communicate stuff when parallel */
if (PAR(cr) && inputrec->eI != eiNM)
global_stat(fplog,gstat,cr,enerd,force_vir,shake_vir,mu_tot,
inputrec,NULL,NULL,NULL,1,&terminate,
top_global,&ems->s,FALSE,
- CGLO_ENERGY |
- CGLO_PRESSURE |
- CGLO_CONSTRAINT |
+ CGLO_ENERGY |
+ CGLO_PRESSURE |
+ CGLO_CONSTRAINT |
CGLO_FIRSTITERATE);
wallcycle_stop(wcycle,ewcMoveE);
}
/* Calculate long range corrections to pressure and energy */
- calc_dispcorr(fplog,inputrec,fr,count,top_global->natoms,ems->s.box,ems->s.lambda,
+ calc_dispcorr(fplog,inputrec,fr,count,top_global->natoms,ems->s.box,ems->s.lambda[efptVDW],
pres,force_vir,&prescorr,&enercorr,&dvdlcorr);
enerd->term[F_DISPCORR] = enercorr;
enerd->term[F_EPOT] += enercorr;
enerd->term[F_DVDL] += dvdlcorr;
ems->epot = enerd->term[F_EPOT];
-
+
if (constr) {
/* Project out the constraint components of the force */
wallcycle_start(wcycle,ewcCONSTR);
- dvdl = 0;
+ dvdlambda = 0;
constrain(NULL,FALSE,FALSE,constr,&top->idef,
inputrec,NULL,cr,count,0,mdatoms,
- ems->s.x,ems->f,ems->f,ems->s.box,ems->s.lambda,&dvdl,
+ ems->s.x,ems->f,ems->f,ems->s.box,ems->s.lambda[efptBONDED],&dvdlambda,
NULL,&shake_vir,nrnb,econqForceDispl,FALSE,0,0);
if (fr->bSepDVDL && fplog)
- fprintf(fplog,sepdvdlformat,"Constraints",t,dvdl);
- enerd->term[F_DHDL_CON] += dvdl;
+ fprintf(fplog,sepdvdlformat,"Constraints",t,dvdlambda);
+ enerd->term[F_DVDL_BONDED] += dvdlambda;
m_add(force_vir,shake_vir,vir);
wallcycle_stop(wcycle,ewcCONSTR);
} else {
clear_mat(ekin);
enerd->term[F_PRES] =
- calc_pres(fr->ePBC,inputrec->nwall,ems->s.box,ekin,vir,pres,
- (fr->eeltype==eelPPPM)?enerd->term[F_COUL_RECIP]:0.0);
+ calc_pres(fr->ePBC,inputrec->nwall,ems->s.box,ekin,vir,pres);
- sum_dhdl(enerd,ems->s.lambda,inputrec);
+ sum_dhdl(enerd,ems->s.lambda,inputrec->fepvals);
if (EI_ENERGY_MINIMIZATION(inputrec->eI))
{
* but to fully optimize this a much more complicated algorithm is required.
*/
snew(fmg,mtop->natoms);
-
+
ncg = s_min->s.ncg_gl;
cg_gl = s_min->s.cg_gl;
i = 0;
i++;
}
}
-
+
sfree(fmg);
return partsum;
* it looks a bit complicated since we take freeze groups into account,
* and might have to sum it in parallel runs.
*/
-
+
if (!DOMAINDECOMP(cr) ||
(s_min->s.ddp_count == cr->dd->ddp_count &&
s_b->s.ddp_count == cr->dd->ddp_count)) {
for(m=0; m<DIM; m++)
if (!opts->nFreeze[gf][m]) {
sum += (fb[i][m] - fm[i][m])*fb[i][m];
- }
+ }
}
} else {
/* We need to reorder cgs while summing */
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,
t_forcerec *fr,
- int repl_ex_nst,int repl_ex_seed,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
rvec *f_global,*p,*sf,*sfm;
double gpa,gpb,gpc,tmp,sum[2],minstep;
real fnormn;
- real stepsize;
+ real stepsize;
real a,b,c,beta=0.0;
real epot_repl=0;
real pnorm;
int number_steps,neval=0,nstcg=inputrec->nstcgsteep;
gmx_mdoutf_t *outf;
int i,m,gf,step,nminstep;
- real terminate=0;
+ real terminate=0;
step=0;
state_global,top_global,s_min,&top,&f,&f_global,
nrnb,mu_tot,fr,&enerd,&graph,mdatoms,&gstat,vsite,constr,
nfile,fnm,&outf,&mdebin);
-
+
/* Print to log file */
print_em_start(fplog,cr,runtime,wcycle,CG);
-
+
/* Max number of steps */
number_steps=inputrec->nsteps;
if (MASTER(cr)) {
/* Copy stuff to the energy bin for easy printing etc. */
upd_mdebin(mdebin,FALSE,FALSE,(double)step,
- mdatoms->tmass,enerd,&s_min->s,s_min->s.box,
- NULL,NULL,vir,pres,NULL,mu_tot,constr);
-
- print_ebin_header(fplog,step,step,s_min->s.lambda);
+ mdatoms->tmass,enerd,&s_min->s,inputrec->fepvals,inputrec->expandedvals,s_min->s.box,
+ NULL,NULL,vir,pres,NULL,mu_tot,constr);
+
+ print_ebin_header(fplog,step,step,s_min->s.lambda[efptFEP]);
print_ebin(outf->fp_ene,TRUE,FALSE,FALSE,fplog,step,step,eprNORMAL,
TRUE,mdebin,fcd,&(top_global->groups),&(inputrec->opts));
}
/* Estimate/guess the initial stepsize */
stepsize = inputrec->em_stepsize/s_min->fnorm;
-
+
if (MASTER(cr)) {
fprintf(stderr," F-max = %12.5e on atom %d\n",
s_min->fmax,s_min->a_fmax+1);
fprintf(fplog," F-Norm = %12.5e\n",
s_min->fnorm/sqrt(state_global->natoms));
fprintf(fplog,"\n");
- }
- /* Start the loop over CG steps.
+ }
+ /* Start the loop over CG steps.
* Each successful step is counted, and we continue until
* we either converge or reach the max number of steps.
*/
converged = FALSE;
for(step=0; (number_steps<0 || (number_steps>=0 && step<=number_steps)) && !converged;step++) {
-
- /* start taking steps in a new direction
- * First time we enter the routine, beta=0, and the direction is
+
+ /* start taking steps in a new direction
+ * First time we enter the routine, beta=0, and the direction is
* simply the negative gradient.
*/
gpa = 0;
gf = 0;
for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++) {
- if (mdatoms->cFREEZE)
+ if (mdatoms->cFREEZE)
gf = mdatoms->cFREEZE[i];
for(m=0; m<DIM; m++) {
if (!inputrec->opts.nFreeze[gf][m]) {
}
}
}
-
+
/* Sum the gradient along the line across CPUs */
if (PAR(cr))
gmx_sumd(1,&gpa,cr);
/* Calculate the norm of the search vector */
get_f_norm_max(cr,&(inputrec->opts),mdatoms,p,&pnorm,NULL,NULL);
-
+
/* Just in case stepsize reaches zero due to numerical precision... */
- if(stepsize<=0)
+ if(stepsize<=0)
stepsize = inputrec->em_stepsize/pnorm;
-
- /*
+
+ /*
* Double check the value of the derivative in the search direction.
* If it is positive it must be due to the old information in the
* CG formula, so just remove that and start over with beta=0.
converged=TRUE;
break;
}
-
+
/* Write coordinates if necessary */
do_x = do_per_step(step,inputrec->nstxout);
do_f = do_per_step(step,inputrec->nstfout);
-
+
write_em_traj(fplog,cr,outf,do_x,do_f,NULL,
top_global,inputrec,step,
s_min,state_global,f_global);
-
+
/* Take a step downhill.
* In theory, we should minimize the function along this direction.
* That is quite possible, but it turns out to take 5-10 function evaluations
* the continue straight to the next CG step without trying to find any minimum.
* If it didn't work (higher energy), there must be a minimum somewhere between
* the old position and the new one.
- *
+ *
* Due to the finite numerical accuracy, it turns out that it is a good idea
* to even accept a SMALL increase in energy, if the derivative is still downhill.
- * This leads to lower final energies in the tests I've done. / Erik
+ * This leads to lower final energies in the tests I've done. / Erik
*/
s_a->epot = s_min->epot;
a = 0.0;
c = a + stepsize; /* reference position along line is zero */
-
+
if (DOMAINDECOMP(cr) && s_min->s.ddp_count < cr->dd->ddp_count) {
em_dd_partition_system(fplog,step,cr,top_global,inputrec,
s_min,top,mdatoms,fr,vsite,constr,
/* Take a trial step (new coords in s_c) */
do_em_step(cr,inputrec,mdatoms,s_min,c,s_min->s.cg_p,s_c,
constr,top,nrnb,wcycle,-1);
-
+
neval++;
/* Calculate energy for the trial step */
evaluate_energy(fplog,bVerbose,cr,
inputrec,nrnb,wcycle,gstat,
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,-1,FALSE);
-
+
/* Calc derivative along line */
p = s_c->s.cg_p;
sf = s_c->f;
gpc=0;
for(i=mdatoms->start; i<mdatoms->start+mdatoms->homenr; i++) {
- for(m=0; m<DIM; m++)
+ for(m=0; m<DIM; m++)
gpc -= p[i][m]*sf[i][m]; /* f is negative gradient, thus the sign */
}
/* Sum the gradient along the line across CPUs */
*/
foundlower = FALSE;
stepsize *= 0.618034;
- }
+ }
+
-
/* OK, if we didn't find a lower value we will have to locate one now - there must
* be one in the interval [a=0,c].
* The same thing is valid here, though: Don't spend dozens of iterations to find
if(gpa<0 && gpc>0)
b = a + gpa*(a-c)/(gpc-gpa);
else
- b = 0.5*(a+c);
-
+ b = 0.5*(a+c);
+
/* safeguard if interpolation close to machine accuracy causes errors:
* never go outside the interval
*/
if(b<=a || b>=c)
b = 0.5*(a+c);
-
+
if (DOMAINDECOMP(cr) && s_min->s.ddp_count != cr->dd->ddp_count) {
/* Reload the old state */
em_dd_partition_system(fplog,-1,cr,top_global,inputrec,
/* Take a trial step to this new point - new coords in s_b */
do_em_step(cr,inputrec,mdatoms,s_min,b,s_min->s.cg_p,s_b,
constr,top,nrnb,wcycle,-1);
-
+
neval++;
/* Calculate energy for the trial step */
evaluate_energy(fplog,bVerbose,cr,
inputrec,nrnb,wcycle,gstat,
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,-1,FALSE);
-
+
/* p does not change within a step, but since the domain decomposition
* might change, we have to use cg_p of s_b here.
*/
/* Sum the gradient along the line across CPUs */
if (PAR(cr))
gmx_sumd(1,&gpb,cr);
-
+
if (debug)
fprintf(debug,"CGE: EpotA %f EpotB %f EpotC %f gpb %f\n",
s_a->epot,s_b->epot,s_c->epot,gpb);
epot_repl = s_b->epot;
-
+
/* Keep one of the intervals based on the value of the derivative at the new point */
if (gpb > 0) {
/* Replace c endpoint with b */
a = b;
gpa = gpb;
}
-
- /*
+
+ /*
* Stop search as soon as we find a value smaller than the endpoints.
* Never run more than 20 steps, no matter what.
*/
nminstep++;
} while ((epot_repl > s_a->epot || epot_repl > s_c->epot) &&
- (nminstep < 20));
-
+ (nminstep < 20));
+
if (fabs(epot_repl - s_min->epot) < fabs(s_min->epot)*GMX_REAL_EPS ||
nminstep >= 20) {
/* OK. We couldn't find a significantly lower energy.
continue;
}
}
-
+
/* Select min energy state of A & C, put the best in B.
*/
if (s_c->epot < s_a->epot) {
gpb = gpa;
b = a;
}
-
+
} else {
if (debug)
fprintf(debug,"CGE: Found a lower energy %f, moving C to B\n",
gpb = gpc;
b = c;
}
-
+
/* new search direction */
/* beta = 0 means forget all memory and restart with steepest descents. */
- if (nstcg && ((step % nstcg)==0))
+ if (nstcg && ((step % nstcg)==0))
beta = 0.0;
else {
/* s_min->fnorm cannot be zero, because then we would have converged
/* Limit beta to prevent oscillations */
if (fabs(beta) > 5.0)
beta = 0.0;
-
-
+
+
/* update positions */
swap_em_state(s_min,s_b);
gpa = gpb;
-
+
/* Print it if necessary */
if (MASTER(cr)) {
if(bVerbose)
s_min->fmax,s_min->a_fmax+1);
/* Store the new (lower) energies */
upd_mdebin(mdebin,FALSE,FALSE,(double)step,
- mdatoms->tmass,enerd,&s_min->s,s_min->s.box,
- NULL,NULL,vir,pres,NULL,mu_tot,constr);
+ mdatoms->tmass,enerd,&s_min->s,inputrec->fepvals,inputrec->expandedvals,s_min->s.box,
+ NULL,NULL,vir,pres,NULL,mu_tot,constr);
+
do_log = do_per_step(step,inputrec->nstlog);
do_ene = do_per_step(step,inputrec->nstenergy);
if(do_log)
- print_ebin_header(fplog,step,step,s_min->s.lambda);
+ print_ebin_header(fplog,step,step,s_min->s.lambda[efptFEP]);
print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,
do_log ? fplog : NULL,step,step,eprNORMAL,
TRUE,mdebin,fcd,&(top_global->groups),&(inputrec->opts));
}
-
+
/* Stop when the maximum force lies below tolerance.
* If we have reached machine precision, converged is already set to true.
- */
+ */
converged = converged || (s_min->fmax < inputrec->em_tol);
-
+
} /* End of the loop */
-
- if (converged)
+
+ if (converged)
step--; /* we never took that last step in this case */
-
+
if (s_min->fmax > inputrec->em_tol)
{
if (MASTER(cr))
warn_step(stderr,inputrec->em_tol,step-1==number_steps,FALSE);
warn_step(fplog ,inputrec->em_tol,step-1==number_steps,FALSE);
}
- converged = FALSE;
+ converged = FALSE;
}
-
+
if (MASTER(cr)) {
/* If we printed energy and/or logfile last step (which was the last step)
* we don't have to do it again, but otherwise print the final values.
*/
if(!do_log) {
/* Write final value to log since we didn't do anything the last step */
- print_ebin_header(fplog,step,step,s_min->s.lambda);
+ print_ebin_header(fplog,step,step,s_min->s.lambda[efptFEP]);
}
if (!do_ene || !do_log) {
/* Write final energy file entries */
/* Print some stuff... */
if (MASTER(cr))
fprintf(stderr,"\nwriting lowest energy coordinates.\n");
-
+
/* IMPORTANT!
* For accurate normal mode calculation it is imperative that we
* store the last conformation into the full precision binary trajectory.
*
* However, we should only do it if we did NOT already write this step
* above (which we did if do_x or do_f was true).
- */
+ */
do_x = !do_per_step(step,inputrec->nstxout);
do_f = (inputrec->nstfout > 0 && !do_per_step(step,inputrec->nstfout));
-
+
write_em_traj(fplog,cr,outf,do_x,do_f,ftp2fn(efSTO,nfile,fnm),
top_global,inputrec,step,
s_min,state_global,f_global);
-
+
fnormn = s_min->fnorm/sqrt(state_global->natoms);
-
+
if (MASTER(cr)) {
print_converged(stderr,CG,inputrec->em_tol,step,converged,number_steps,
s_min->epot,s_min->fmax,s_min->a_fmax,fnormn);
print_converged(fplog,CG,inputrec->em_tol,step,converged,number_steps,
s_min->epot,s_min->fmax,s_min->a_fmax,fnormn);
-
+
fprintf(fplog,"\nPerformed %d energy evaluations in total.\n",neval);
}
-
+
finish_em(fplog,cr,outf,runtime,wcycle);
-
+
/* To print the actual number of steps we needed somewhere */
runtime->nsteps_done = step;
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,
t_forcerec *fr,
- int repl_ex_nst,int repl_ex_seed,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
rvec *f_global;
int ncorr,nmaxcorr,point,cp,neval,nminstep;
double stepsize,gpa,gpb,gpc,tmp,minstep;
- real *rho,*alpha,*ff,*xx,*p,*s,*lastx,*lastf,**dx,**dg;
+ real *rho,*alpha,*ff,*xx,*p,*s,*lastx,*lastf,**dx,**dg;
real *xa,*xb,*xc,*fa,*fb,*fc,*xtmp,*ftmp;
real a,b,c,maxdelta,delta;
real diag,Epot0,Epot,EpotA,EpotB,EpotC;
int start,end,number_steps;
gmx_mdoutf_t *outf;
int i,k,m,n,nfmax,gf,step;
+ int mdof_flags;
/* not used */
real terminate;
if (PAR(cr))
gmx_fatal(FARGS,"Cannot do parallel L-BFGS Minimization - yet.\n");
-
+
n = 3*state->natoms;
nmaxcorr = inputrec->nbfgscorr;
-
+
/* Allocate memory */
/* Use pointers to real so we dont have to loop over both atoms and
* dimensions all the time...
snew(fc,n);
snew(frozen,n);
- snew(p,n);
- snew(lastx,n);
- snew(lastf,n);
+ snew(p,n);
+ snew(lastx,n);
+ snew(lastf,n);
snew(rho,nmaxcorr);
snew(alpha,nmaxcorr);
-
+
snew(dx,nmaxcorr);
for(i=0;i<nmaxcorr;i++)
snew(dx[i],n);
-
+
snew(dg,nmaxcorr);
for(i=0;i<nmaxcorr;i++)
snew(dg[i],n);
step = 0;
- neval = 0;
+ neval = 0;
/* Init em */
init_em(fplog,LBFGS,cr,inputrec,
start = mdatoms->start;
end = mdatoms->homenr + start;
-
+
/* Print to log file */
print_em_start(fplog,cr,runtime,wcycle,LBFGS);
-
+
do_log = do_ene = do_x = do_f = TRUE;
-
+
/* Max number of steps */
number_steps=inputrec->nsteps;
for(i=start; i<end; i++) {
if (mdatoms->cFREEZE)
gf = mdatoms->cFREEZE[i];
- for(m=0; m<DIM; m++)
- frozen[3*i+m]=inputrec->opts.nFreeze[gf][m];
+ for(m=0; m<DIM; m++)
+ frozen[3*i+m]=inputrec->opts.nFreeze[gf][m];
}
if (MASTER(cr))
sp_header(stderr,LBFGS,inputrec->em_tol,number_steps);
if (fplog)
sp_header(fplog,LBFGS,inputrec->em_tol,number_steps);
-
+
if (vsite)
construct_vsites(fplog,vsite,state->x,nrnb,1,NULL,
top->idef.iparams,top->idef.il,
fr->ePBC,fr->bMolPBC,graph,cr,state->box);
-
+
/* Call the force routine and some auxiliary (neighboursearching etc.) */
/* do_force always puts the charge groups in the box and shifts again
* We do not unshift, so molecules are always whole
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,-1,TRUE);
where();
-
+
if (MASTER(cr)) {
/* Copy stuff to the energy bin for easy printing etc. */
upd_mdebin(mdebin,FALSE,FALSE,(double)step,
- mdatoms->tmass,enerd,state,state->box,
- NULL,NULL,vir,pres,NULL,mu_tot,constr);
-
- print_ebin_header(fplog,step,step,state->lambda);
+ mdatoms->tmass,enerd,state,inputrec->fepvals,inputrec->expandedvals,state->box,
+ NULL,NULL,vir,pres,NULL,mu_tot,constr);
+
+ print_ebin_header(fplog,step,step,state->lambda[efptFEP]);
print_ebin(outf->fp_ene,TRUE,FALSE,FALSE,fplog,step,step,eprNORMAL,
TRUE,mdebin,fcd,&(top_global->groups),&(inputrec->opts));
}
where();
-
+
/* This is the starting energy */
Epot = enerd->term[F_EPOT];
-
+
fnorm = ems.fnorm;
fmax = ems.fmax;
nfmax = ems.a_fmax;
-
+
/* Set the initial step.
- * since it will be multiplied by the non-normalized search direction
+ * since it will be multiplied by the non-normalized search direction
* vector (force vector the first time), we scale it by the
* norm of the force.
*/
-
+
if (MASTER(cr)) {
fprintf(stderr,"Using %d BFGS correction steps.\n\n",nmaxcorr);
fprintf(stderr," F-max = %12.5e on atom %d\n",fmax,nfmax+1);
fprintf(fplog," F-max = %12.5e on atom %d\n",fmax,nfmax+1);
fprintf(fplog," F-Norm = %12.5e\n",fnorm/sqrt(state->natoms));
fprintf(fplog,"\n");
- }
-
+ }
+
point=0;
for(i=0;i<n;i++)
if(!frozen[i])
stepsize = 1.0/fnorm;
converged = FALSE;
-
- /* Start the loop over BFGS steps.
+
+ /* Start the loop over BFGS steps.
* Each successful step is counted, and we continue until
* we either converge or reach the max number of steps.
*/
-
+
ncorr=0;
/* Set the gradient from the force */
converged = FALSE;
for(step=0; (number_steps<0 || (number_steps>=0 && step<=number_steps)) && !converged; step++) {
-
+
/* Write coordinates if necessary */
do_x = do_per_step(step,inputrec->nstxout);
do_f = do_per_step(step,inputrec->nstfout);
-
+
- write_traj(fplog,cr,outf,MDOF_X | MDOF_F,
+ mdof_flags = 0;
+ if (do_x)
+ {
+ mdof_flags |= MDOF_X;
+ }
+
+ if (do_f)
+ {
+ mdof_flags |= MDOF_F;
+ }
+
+ write_traj(fplog,cr,outf,mdof_flags,
top_global,step,(real)step,state,state,f,f,NULL,NULL);
/* Do the linesearching in the direction dx[point][0..(n-1)] */
-
+
/* pointer to current direction - point=0 first time here */
s=dx[point];
-
+
/* calculate line gradient */
- for(gpa=0,i=0;i<n;i++)
+ for(gpa=0,i=0;i<n;i++)
gpa-=s[i]*ff[i];
- /* Calculate minimum allowed stepsize, before the average (norm)
- * relative change in coordinate is smaller than precision
+ /* Calculate minimum allowed stepsize, before the average (norm)
+ * relative change in coordinate is smaller than precision
*/
for(minstep=0,i=0;i<n;i++) {
tmp=fabs(xx[i]);
minstep += tmp*tmp;
}
minstep = GMX_REAL_EPS/sqrt(minstep/n);
-
+
if(stepsize<minstep) {
converged=TRUE;
break;
}
-
+
/* Store old forces and coordinates */
for(i=0;i<n;i++) {
lastx[i]=xx[i];
lastf[i]=ff[i];
}
Epot0=Epot;
-
+
first=TRUE;
-
+
for(i=0;i<n;i++)
xa[i]=xx[i];
-
+
/* Take a step downhill.
* In theory, we should minimize the function along this direction.
* That is quite possible, but it turns out to take 5-10 function evaluations
* the continue straight to the next BFGS step without trying to find any minimum.
* If it didn't work (higher energy), there must be a minimum somewhere between
* the old position and the new one.
- *
+ *
* Due to the finite numerical accuracy, it turns out that it is a good idea
* to even accept a SMALL increase in energy, if the derivative is still downhill.
- * This leads to lower final energies in the tests I've done. / Erik
+ * This leads to lower final energies in the tests I've done. / Erik
*/
foundlower=FALSE;
EpotA = Epot0;
a = 0.0;
c = a + stepsize; /* reference position along line is zero */
- /* Check stepsize first. We do not allow displacements
+ /* Check stepsize first. We do not allow displacements
* larger than emstep.
*/
do {
/* Take a trial step */
for (i=0; i<n; i++)
xc[i] = lastx[i] + c*s[i];
-
+
neval++;
/* Calculate energy for the trial step */
ems.s.x = (rvec *)xc;
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,step,FALSE);
EpotC = ems.epot;
-
+
/* Calc derivative along line */
for(gpc=0,i=0; i<n; i++) {
gpc -= s[i]*fc[i]; /* f is negative gradient, thus the sign */
/* Sum the gradient along the line across CPUs */
if (PAR(cr))
gmx_sumd(1,&gpc,cr);
-
+
/* This is the max amount of increase in energy we tolerate */
tmp=sqrt(GMX_REAL_EPS)*fabs(EpotA);
-
+
/* Accept the step if the energy is lower, or if it is not significantly higher
* and the line derivative is still negative.
*/
*/
foundlower = FALSE;
stepsize *= 0.618034;
- }
-
+ }
+
/* OK, if we didn't find a lower value we will have to locate one now - there must
* be one in the interval [a=0,c].
* The same thing is valid here, though: Don't spend dozens of iterations to find
*/
if(!foundlower) {
-
+
nminstep=0;
do {
/* Select a new trial point.
* If the derivatives at points a & c have different sign we interpolate to zero,
* otherwise just do a bisection.
*/
-
+
if(gpa<0 && gpc>0)
b = a + gpa*(a-c)/(gpc-gpa);
else
- b = 0.5*(a+c);
-
+ b = 0.5*(a+c);
+
/* safeguard if interpolation close to machine accuracy causes errors:
* never go outside the interval
*/
if(b<=a || b>=c)
b = 0.5*(a+c);
-
+
/* Take a trial step */
- for (i=0; i<n; i++)
+ for (i=0; i<n; i++)
xb[i] = lastx[i] + b*s[i];
-
+
neval++;
/* Calculate energy for the trial step */
ems.s.x = (rvec *)xb;
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,step,FALSE);
EpotB = ems.epot;
-
+
fnorm = ems.fnorm;
-
- for(gpb=0,i=0; i<n; i++)
+
+ for(gpb=0,i=0; i<n; i++)
gpb -= s[i]*fb[i]; /* f is negative gradient, thus the sign */
-
+
/* Sum the gradient along the line across CPUs */
if (PAR(cr))
gmx_sumd(1,&gpb,cr);
-
+
/* Keep one of the intervals based on the value of the derivative at the new point */
if(gpb>0) {
/* Replace c endpoint with b */
c = b;
gpc = gpb;
/* swap coord pointers b/c */
- xtmp = xb;
+ xtmp = xb;
ftmp = fb;
- xb = xc;
+ xb = xc;
fb = fc;
xc = xtmp;
fc = ftmp;
a = b;
gpa = gpb;
/* swap coord pointers a/b */
- xtmp = xb;
+ xtmp = xb;
ftmp = fb;
- xb = xa;
+ xb = xa;
fb = fa;
- xa = xtmp;
+ xa = xtmp;
fa = ftmp;
}
-
- /*
+
+ /*
* Stop search as soon as we find a value smaller than the endpoints,
* or if the tolerance is below machine precision.
* Never run more than 20 steps, no matter what.
*/
- nminstep++;
+ nminstep++;
} while((EpotB>EpotA || EpotB>EpotC) && (nminstep<20));
if(fabs(EpotB-Epot0)<GMX_REAL_EPS || nminstep>=20) {
continue;
}
}
-
+
/* Select min energy state of A & C, put the best in xx/ff/Epot
*/
if(EpotC<EpotA) {
}
stepsize=a;
}
-
+
} else {
/* found lower */
Epot = EpotC;
stepsize=c;
}
- /* Update the memory information, and calculate a new
- * approximation of the inverse hessian
+ /* Update the memory information, and calculate a new
+ * approximation of the inverse hessian
*/
-
- /* Have new data in Epot, xx, ff */
+
+ /* Have new data in Epot, xx, ff */
if(ncorr<nmaxcorr)
ncorr++;
dg[point][i]=lastf[i]-ff[i];
dx[point][i]*=stepsize;
}
-
+
dgdg=0;
- dgdx=0;
+ dgdx=0;
for(i=0;i<n;i++) {
dgdg+=dg[point][i]*dg[point][i];
dgdx+=dg[point][i]*dx[point][i];
}
-
+
diag=dgdx/dgdg;
-
+
rho[point]=1.0/dgdx;
point++;
-
+
if(point>=nmaxcorr)
point=0;
-
+
/* Update */
for(i=0;i<n;i++)
p[i]=ff[i];
-
+
cp=point;
-
+
/* Recursive update. First go back over the memory points */
for(k=0;k<ncorr;k++) {
cp--;
- if(cp<0)
+ if(cp<0)
cp=ncorr-1;
-
+
sq=0;
for(i=0;i<n;i++)
sq+=dx[cp][i]*p[i];
-
+
alpha[cp]=rho[cp]*sq;
-
+
for(i=0;i<n;i++)
- p[i] -= alpha[cp]*dg[cp][i];
+ p[i] -= alpha[cp]*dg[cp][i];
}
-
+
for(i=0;i<n;i++)
p[i] *= diag;
-
+
/* And then go forward again */
for(k=0;k<ncorr;k++) {
yr = 0;
for(i=0;i<n;i++)
yr += p[i]*dg[cp][i];
-
- beta = rho[cp]*yr;
+
+ beta = rho[cp]*yr;
beta = alpha[cp]-beta;
-
+
for(i=0;i<n;i++)
p[i] += beta*dx[cp][i];
-
- cp++;
+
+ cp++;
if(cp>=ncorr)
cp=0;
}
-
+
for(i=0;i<n;i++)
if(!frozen[i])
dx[point][i] = p[i];
dx[point][i] = 0;
stepsize=1.0;
-
+
/* Test whether the convergence criterion is met */
get_f_norm_max(cr,&(inputrec->opts),mdatoms,f,&fnorm,&fmax,&nfmax);
-
+
/* Print it if necessary */
if (MASTER(cr)) {
if(bVerbose)
step,Epot,fnorm/sqrt(state->natoms),fmax,nfmax+1);
/* Store the new (lower) energies */
upd_mdebin(mdebin,FALSE,FALSE,(double)step,
- mdatoms->tmass,enerd,state,state->box,
- NULL,NULL,vir,pres,NULL,mu_tot,constr);
+ mdatoms->tmass,enerd,state,inputrec->fepvals,inputrec->expandedvals,state->box,
+ NULL,NULL,vir,pres,NULL,mu_tot,constr);
do_log = do_per_step(step,inputrec->nstlog);
do_ene = do_per_step(step,inputrec->nstenergy);
if(do_log)
- print_ebin_header(fplog,step,step,state->lambda);
+ print_ebin_header(fplog,step,step,state->lambda[efptFEP]);
print_ebin(outf->fp_ene,do_ene,FALSE,FALSE,
do_log ? fplog : NULL,step,step,eprNORMAL,
TRUE,mdebin,fcd,&(top_global->groups),&(inputrec->opts));
}
-
+
/* Stop when the maximum force lies below tolerance.
* If we have reached machine precision, converged is already set to true.
*/
-
+
converged = converged || (fmax < inputrec->em_tol);
-
+
} /* End of the loop */
-
- if(converged)
+
+ if(converged)
step--; /* we never took that last step in this case */
-
+
if(fmax>inputrec->em_tol)
{
if (MASTER(cr))
warn_step(stderr,inputrec->em_tol,step-1==number_steps,FALSE);
warn_step(fplog ,inputrec->em_tol,step-1==number_steps,FALSE);
}
- converged = FALSE;
+ converged = FALSE;
}
-
+
/* If we printed energy and/or logfile last step (which was the last step)
* we don't have to do it again, but otherwise print the final values.
*/
if(!do_log) /* Write final value to log since we didn't do anythin last step */
- print_ebin_header(fplog,step,step,state->lambda);
+ print_ebin_header(fplog,step,step,state->lambda[efptFEP]);
if(!do_ene || !do_log) /* Write final energy file entries */
print_ebin(outf->fp_ene,!do_ene,FALSE,FALSE,
!do_log ? fplog : NULL,step,step,eprNORMAL,
TRUE,mdebin,fcd,&(top_global->groups),&(inputrec->opts));
-
+
/* Print some stuff... */
if (MASTER(cr))
fprintf(stderr,"\nwriting lowest energy coordinates.\n");
-
+
/* IMPORTANT!
* For accurate normal mode calculation it is imperative that we
* store the last conformation into the full precision binary trajectory.
*
* However, we should only do it if we did NOT already write this step
* above (which we did if do_x or do_f was true).
- */
+ */
do_x = !do_per_step(step,inputrec->nstxout);
do_f = !do_per_step(step,inputrec->nstfout);
write_em_traj(fplog,cr,outf,do_x,do_f,ftp2fn(efSTO,nfile,fnm),
top_global,inputrec,step,
&ems,state,f);
-
+
if (MASTER(cr)) {
print_converged(stderr,LBFGS,inputrec->em_tol,step,converged,
number_steps,Epot,fmax,nfmax,fnorm/sqrt(state->natoms));
print_converged(fplog,LBFGS,inputrec->em_tol,step,converged,
number_steps,Epot,fmax,nfmax,fnorm/sqrt(state->natoms));
-
+
fprintf(fplog,"\nPerformed %d energy evaluations in total.\n",neval);
}
-
+
finish_em(fplog,cr,outf,runtime,wcycle);
/* To print the actual number of steps we needed somewhere */
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,
t_forcerec *fr,
- int repl_ex_nst,int repl_ex_seed,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
gmx_runtime_t *runtime)
-{
+{
const char *SD="Steepest Descents";
em_state_t *s_min,*s_try;
rvec *f_global;
real stepsize,constepsize;
real ustep,dvdlambda,fnormn;
gmx_mdoutf_t *outf;
- t_mdebin *mdebin;
- gmx_bool bDone,bAbort,do_x,do_f;
- tensor vir,pres;
+ t_mdebin *mdebin;
+ gmx_bool bDone,bAbort,do_x,do_f;
+ tensor vir,pres;
rvec mu_tot;
int nsteps;
- int count=0;
- int steps_accepted=0;
+ int count=0;
+ int steps_accepted=0;
/* not used */
real terminate=0;
state_global,top_global,s_try,&top,&f,&f_global,
nrnb,mu_tot,fr,&enerd,&graph,mdatoms,&gstat,vsite,constr,
nfile,fnm,&outf,&mdebin);
-
+
/* Print to log file */
print_em_start(fplog,cr,runtime,wcycle,SD);
-
- /* Set variables for stepsize (in nm). This is the largest
- * step that we are going to make in any direction.
+
+ /* Set variables for stepsize (in nm). This is the largest
+ * step that we are going to make in any direction.
*/
- ustep = inputrec->em_stepsize;
+ ustep = inputrec->em_stepsize;
stepsize = 0;
-
+
/* Max number of steps */
- nsteps = inputrec->nsteps;
-
- if (MASTER(cr))
+ nsteps = inputrec->nsteps;
+
+ if (MASTER(cr))
/* Print to the screen */
sp_header(stderr,SD,inputrec->em_tol,nsteps);
if (fplog)
sp_header(fplog,SD,inputrec->em_tol,nsteps);
-
+
/**** HERE STARTS THE LOOP ****
- * count is the counter for the number of steps
+ * count is the counter for the number of steps
* bDone will be TRUE when the minimization has converged
* bAbort will be TRUE when nsteps steps have been performed or when
* the stepsize becomes smaller than is reasonable for machine precision
bAbort = FALSE;
while( !bDone && !bAbort ) {
bAbort = (nsteps >= 0) && (count == nsteps);
-
+
/* set new coordinates, except for first step */
if (count > 0) {
do_em_step(cr,inputrec,mdatoms,s_min,stepsize,s_min->f,s_try,
constr,top,nrnb,wcycle,count);
}
-
+
evaluate_energy(fplog,bVerbose,cr,
state_global,top_global,s_try,top,
inputrec,nrnb,wcycle,gstat,
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,count,count==0);
-
+
if (MASTER(cr))
- print_ebin_header(fplog,count,count,s_try->s.lambda);
+ print_ebin_header(fplog,count,count,s_try->s.lambda[efptFEP]);
if (count == 0)
s_min->epot = s_try->epot + 1;
-
+
/* Print it if necessary */
if (MASTER(cr)) {
if (bVerbose) {
count,ustep,s_try->epot,s_try->fmax,s_try->a_fmax+1,
(s_try->epot < s_min->epot) ? '\n' : '\r');
}
-
+
if (s_try->epot < s_min->epot) {
/* Store the new (lower) energies */
upd_mdebin(mdebin,FALSE,FALSE,(double)count,
- mdatoms->tmass,enerd,&s_try->s,s_try->s.box,
- NULL,NULL,vir,pres,NULL,mu_tot,constr);
+ mdatoms->tmass,enerd,&s_try->s,inputrec->fepvals,inputrec->expandedvals,
+ s_try->s.box, NULL,NULL,vir,pres,NULL,mu_tot,constr);
print_ebin(outf->fp_ene,TRUE,
do_per_step(steps_accepted,inputrec->nstdisreout),
do_per_step(steps_accepted,inputrec->nstorireout),
mdebin,fcd,&(top_global->groups),&(inputrec->opts));
fflush(fplog);
}
- }
-
- /* Now if the new energy is smaller than the previous...
+ }
+
+ /* Now if the new energy is smaller than the previous...
* or if this is the first step!
- * or if we did random steps!
+ * or if we did random steps!
*/
-
+
if ( (count==0) || (s_try->epot < s_min->epot) ) {
- steps_accepted++;
+ steps_accepted++;
/* Test whether the convergence criterion is met... */
bDone = (s_try->fmax < inputrec->em_tol);
-
+
/* Copy the arrays for force, positions and energy */
- /* The 'Min' array always holds the coords and forces of the minimal
+ /* The 'Min' array always holds the coords and forces of the minimal
sampled energy */
swap_em_state(s_min,s_try);
if (count > 0)
write_em_traj(fplog,cr,outf,do_x,do_f,NULL,
top_global,inputrec,count,
s_min,state_global,f_global);
- }
+ }
else {
/* If energy is not smaller make the step smaller... */
ustep *= 0.5;
nrnb,wcycle);
}
}
-
+
/* Determine new step */
stepsize = ustep/s_min->fmax;
-
+
/* Check if stepsize is too small, with 1 nm as a characteristic length */
#ifdef GMX_DOUBLE
if (count == nsteps || ustep < 1e-12)
}
bAbort=TRUE;
}
-
+
count++;
} /* End of the loop */
-
+
/* Print some shit... */
- if (MASTER(cr))
- fprintf(stderr,"\nwriting lowest energy coordinates.\n");
+ if (MASTER(cr))
+ fprintf(stderr,"\nwriting lowest energy coordinates.\n");
write_em_traj(fplog,cr,outf,TRUE,inputrec->nstfout,ftp2fn(efSTO,nfile,fnm),
top_global,inputrec,count,
s_min,state_global,f_global);
}
finish_em(fplog,cr,outf,runtime,wcycle);
-
+
/* To print the actual number of steps we needed somewhere */
inputrec->nsteps=count;
runtime->nsteps_done = count;
-
+
return 0;
} /* That's all folks */
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,
t_forcerec *fr,
- int repl_ex_nst,int repl_ex_seed,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ gmx_membed_t membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
- real t,lambda;
+ real t,t0,lambda,lam0;
gmx_bool bNS;
tensor vir,pres;
rvec mu_tot;
gmx_sparsematrix_t * sparse_matrix = NULL;
real * full_matrix = NULL;
em_state_t * state_work;
-
+
/* added with respect to mdrun */
int i,j,k,row,col;
real der_range=10.0*sqrt(GMX_REAL_EPS);
real x_min;
real fnorm,fmax;
-
+
if (constr != NULL)
{
gmx_fatal(FARGS,"Constraints present with Normal Mode Analysis, this combination is not supported");
}
state_work = init_em_state();
-
+
/* Init em and store the local state in state_minimum */
init_em(fplog,NM,cr,inputrec,
state_global,top_global,state_work,&top,
&f,&f_global,
nrnb,mu_tot,fr,&enerd,&graph,mdatoms,&gstat,vsite,constr,
nfile,fnm,&outf,NULL);
-
+
natoms = top_global->natoms;
snew(fneg,natoms);
snew(dfdx,natoms);
-
+
#ifndef GMX_DOUBLE
if (MASTER(cr))
{
" are fairly modest even if you recompile in double precision.\n\n");
}
#endif
-
+
/* Check if we can/should use sparse storage format.
*
* Sparse format is only useful when the Hessian itself is sparse, which it
- * will be when we use a cutoff.
+ * will be when we use a cutoff.
* For small systems (n<1000) it is easier to always use full matrix format, though.
*/
if(EEL_FULL(fr->eeltype) || fr->rlist==0.0)
fprintf(stderr,"Using compressed symmetric sparse Hessian format.\n");
bSparse = TRUE;
}
-
+
sz = DIM*top_global->natoms;
-
+
fprintf(stderr,"Allocating Hessian memory...\n\n");
if(bSparse)
{
snew(full_matrix,sz*sz);
}
-
+
/* Initial values */
- t = inputrec->init_t;
- lambda = inputrec->init_lambda;
-
+ t0 = inputrec->init_t;
+ lam0 = inputrec->fepvals->init_lambda;
+ t = t0;
+ lambda = lam0;
+
init_nrnb(nrnb);
-
+
where();
-
+
/* Write start time and temperature */
print_em_start(fplog,cr,runtime,wcycle,NM);
/* fudge nr of steps to nr of atoms */
inputrec->nsteps = natoms*2;
- if (MASTER(cr))
+ if (MASTER(cr))
{
fprintf(stderr,"starting normal mode calculation '%s'\n%d steps.\n\n",
*(top_global->name),(int)inputrec->nsteps);
}
nnodes = cr->nnodes;
-
+
/* Make evaluate_energy do a single node force calculation */
cr->nnodes = 1;
evaluate_energy(fplog,bVerbose,cr,
if (MASTER(cr))
{
fprintf(stderr,"Maximum force:%12.5e\n",state_work->fmax);
- if (state_work->fmax > 1.0e-3)
+ if (state_work->fmax > 1.0e-3)
{
fprintf(stderr,"Maximum force probably not small enough to");
fprintf(stderr," ensure that you are in an \nenergy well. ");
fprintf(stderr," when the\nresulting matrix is diagonalized.\n");
}
}
-
+
/***********************************************************
*
- * Loop over all pairs in matrix
- *
- * do_force called twice. Once with positive and
- * once with negative displacement
+ * Loop over all pairs in matrix
+ *
+ * do_force called twice. Once with positive and
+ * once with negative displacement
*
************************************************************/
/* Steps are divided one by one over the nodes */
- for(atom=cr->nodeid; atom<natoms; atom+=nnodes)
+ for(atom=cr->nodeid; atom<natoms; atom+=nnodes)
{
-
- for (d=0; d<DIM; d++)
+
+ for (d=0; d<DIM; d++)
{
x_min = state_work->s.x[atom][d];
state_work->s.x[atom][d] = x_min - der_range;
-
+
/* Make evaluate_energy do a single node force calculation */
cr->nnodes = 1;
evaluate_energy(fplog,bVerbose,cr,
inputrec,nrnb,wcycle,gstat,
vsite,constr,fcd,graph,mdatoms,fr,
mu_tot,enerd,vir,pres,atom*2,FALSE);
-
+
for(i=0; i<natoms; i++)
{
copy_rvec(state_work->f[i], fneg[i]);
}
-
+
state_work->s.x[atom][d] = x_min + der_range;
-
+
evaluate_energy(fplog,bVerbose,cr,
state_global,top_global,state_work,top,
inputrec,nrnb,wcycle,gstat,
/* x is restored to original */
state_work->s.x[atom][d] = x_min;
- for(j=0; j<natoms; j++)
+ for(j=0; j<natoms; j++)
{
- for (k=0; (k<DIM); k++)
+ for (k=0; (k<DIM); k++)
{
dfdx[j][k] =
-(state_work->f[j][k] - fneg[j][k])/(2*der_range);
row = (atom + node)*DIM + d;
- for(j=0; j<natoms; j++)
+ for(j=0; j<natoms; j++)
{
- for(k=0; k<DIM; k++)
+ for(k=0; k<DIM; k++)
{
col = j*DIM + k;
-
+
if (bSparse)
{
if (col >= row && dfdx[j][k] != 0.0)
}
}
}
-
+
if (bVerbose && fplog)
{
- fflush(fplog);
+ fflush(fplog);
}
}
/* write progress */
- if (MASTER(cr) && bVerbose)
+ if (MASTER(cr) && bVerbose)
{
fprintf(stderr,"\rFinished step %d out of %d",
- min(atom+nnodes,natoms),natoms);
+ min(atom+nnodes,natoms),natoms);
fflush(stderr);
}
}
-
- if (MASTER(cr))
+
+ if (MASTER(cr))
{
fprintf(stderr,"\n\nWriting Hessian...\n");
gmx_mtxio_write(ftp2fn(efMTX,nfile,fnm),sz,sz,full_matrix,sparse_matrix);
finish_em(fplog,cr,outf,runtime,wcycle);
runtime->nsteps_done = natoms*2;
-
+
return 0;
}
+
#include <math.h>
/*#define HAVE_NN_LOOPS*/
-/* Set environment variable CFLAGS = "-fopenmp" when running
- * configure and define DOUSEOPENMP to make use of parallelized
- * calculation of autocorrelation function.
- * It also adds a new option -nthreads which sets the number of threads.
- * */
-/*#define DOUSEOPENMP*/
-
-#ifdef DOUSEOPENMP
-#define HAVE_OPENMP
-#include "omp.h"
-#endif
+
+#include "gmx_omp.h"
#include "statutil.h"
#include "copyrite.h"
{"n-n","n-n+1","n-n+2","n-n+3","n-n+4","n-n+5","n-n>6"};
#define MAXHH 4
-#ifdef HAVE_OPENMP
+#ifdef GMX_OPENMP
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==0)
#else
#define MASTER_THREAD_ONLY(threadNr) ((threadNr)==(threadNr))
/* Not found apparently. Add it to the list! */
/* printf("New shift found: %i,%i,%i\n",r[XX],r[YY],r[ZZ]); */
-/* Unfortunately this needs to be critical it seems. */
-#ifdef HAVE_OPENMP
#pragma omp critical
-#endif
{
if (!per->p2i) {
fprintf(stderr, "p2i not initialized. This shouldn't happen!\n");
E = 0;
hb->hbE.E[d][a][h][frame] = E;
-#ifdef HAVE_OPENMP
+
#pragma omp critical
-#endif
{
hb->hbE.Etot[frame] += E;
}
}
}
+
}
static void inc_nhbonds(t_donors *ddd,int d, int h)
grpa,hb->a.grp[ia],a+1);
if (bMerge)
- if ((daSwap = isInterchangable(hb, d, a, grpd, grpa) || bContact) && d>a)
+ {
+
+ if (isInterchangable(hb, d, a, grpd, grpa) && d>a)
/* Then swap identity so that the id of d is lower then that of a.
*
* This should really be redundant by now, as is_hbond() now ought to return
* hbNo in the cases where this conditional is TRUE. */
{
+ daSwap = TRUE;
k = d;
d = a;
a = k;
gmx_fatal(FARGS,"Inconsistent acceptor groups, %d iso %d, atom %d",
grpa,hb->a.grp[ia],a+1);
}
+ }
if (hb->hbmap) {
/* Loop over hydrogens to find which hydrogen is in this particular HB */
k = 0;
if (hb->bHBmap) {
- if (hb->hbmap[id][ia] == NULL) {
- snew(hb->hbmap[id][ia],1);
- snew(hb->hbmap[id][ia]->h,hb->maxhydro);
- snew(hb->hbmap[id][ia]->g,hb->maxhydro);
+
+#pragma omp critical
+ {
+ if (hb->hbmap[id][ia] == NULL) {
+ snew(hb->hbmap[id][ia],1);
+ snew(hb->hbmap[id][ia]->h,hb->maxhydro);
+ snew(hb->hbmap[id][ia]->g,hb->maxhydro);
+ }
+ add_ff(hb,id,k,ia,frame,ihb,p);
}
- add_ff(hb,id,k,ia,frame,ihb,p);
}
/* Strange construction with frame >=0 is a relic from old code
inc_nhbonds(&(hb->d),d,h);
}
-/* Now a redundant function. It might find use at some point though. */
-static gmx_bool in_list(atom_id selection,int isize,atom_id *index)
-{
- int i;
- gmx_bool bFound;
-
- bFound=FALSE;
- for(i=0; (i<isize) && !bFound; i++)
- if(selection == index[i])
- bFound=TRUE;
-
- return bFound;
-}
-
static char *mkatomname(t_atoms *atoms,int i)
{
static char buf[32];
int i,j,nra,n;
t_functype func_type;
t_ilist *interaction;
- atom_id nr1,nr2;
+ atom_id nr1,nr2,nr3;
gmx_bool stop;
if (!ddd->dptr) {
/* check out this functype */
if (func_type == F_SETTLE) {
- nr1=interaction->iatoms[i+1];
+ nr1 = interaction->iatoms[i+1];
+ nr2 = interaction->iatoms[i+2];
+ nr3 = interaction->iatoms[i+3];
if (ISINGRP(datable[nr1])) {
- if (ISINGRP(datable[nr1+1])) {
+ if (ISINGRP(datable[nr2])) {
datable[nr1] |= DON;
add_dh(ddd,nr1,nr1+1,grp,datable);
}
- if (ISINGRP(datable[nr1+2])) {
+ if (ISINGRP(datable[nr3])) {
datable[nr1] |= DON;
add_dh(ddd,nr1,nr1+2,grp,datable);
}
return grid;
}
-static void control_pHist(t_hbdata *hb, int nframes)
-{
- int i,j,k;
- PSTYPE p;
- for (i=0;i<hb->d.nrd;i++)
- for (j=0;j<hb->a.nra;j++)
- if (hb->per->pHist[i][j].len != 0)
- for (k=hb->hbmap[i][j][0].n0; k<nframes; k++) {
- p = getPshift(hb->per->pHist[i][j], k);
- if (p>hb->per->nper)
- fprintf(stderr, "Weird stuff in pHist[%i][%i].p at frame %i: p=%i\n",
- i,j,k,p);
- }
-}
-
static void reset_nhbonds(t_donors *ddd)
{
int i,j;
* This could be implemented slightly more efficient, but the code
* would get much more complicated.
*/
-#define B(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? 0 : (x-1))
-#define E(n,x,bTric,bEdge) ((n==1) ? x : bTric&&(bEdge) ? n-1 : (x+1))
-#define GRIDMOD(j,n) (j+n)%(n)
-#define LOOPGRIDINNER(x,y,z,xx,yy,zz,xo,yo,zo,n,bTric) \
- for(zz=B(n[ZZ],zo,bTric,FALSE); zz<=E(n[ZZ],zo,bTric,FALSE); zz++) { \
- z=GRIDMOD(zz,n[ZZ]); \
- for(yy=B(n[YY],yo,bTric,z==0||z==n[ZZ]-1); \
- yy<=E(n[YY],yo,bTric,z==0||z==n[ZZ]-1); yy++) { \
- y=GRIDMOD(yy,n[YY]); \
- for(xx=B(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); \
- xx<=E(n[XX],xo,bTric,y==0||y==n[YY]-1||z==0||z==n[ZZ]-1); xx++) { \
- x=GRIDMOD(xx,n[XX]);
-#define ENDLOOPGRIDINNER \
- } \
- } \
- } \
- \
+static inline gmx_bool grid_loop_begin(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+{
+ return ((n==1) ? x : bTric && bEdge ? 0 : (x-1));
+}
+static inline gmx_bool grid_loop_end(int n, int x, gmx_bool bTric, gmx_bool bEdge)
+{
+ return ((n==1) ? x : bTric && bEdge ? (n-1) : (x+1));
+}
+static inline int grid_mod(int j, int n)
+{
+ return (j+n) % (n);
+}
static void dump_grid(FILE *fp, ivec ngrid, t_gridcell ***grid)
{
g=NULL;
}
-static void pbc_correct(rvec dx,matrix box,rvec hbox)
-{
- int m;
- for(m=DIM-1; m>=0; m--) {
- if ( dx[m] < -hbox[m] )
- rvec_inc(dx,box[m]);
- else if ( dx[m] >= hbox[m] )
- rvec_dec(dx,box[m]);
- }
-}
-
void pbc_correct_gem(rvec dx,matrix box,rvec hbox)
{
int m;
rvec_sub(x[d],x[a],r_da);
/* Insert projection code here */
- /* if (d>a && ((isInterchangable(hb, d, a, grpd, grpa) && bMerge) || bContact)) */
- /* /\* Then this hbond will be found again, or it has already been found. *\/ */
- /* return hbNo; */
-
+ if (bMerge && d>a && isInterchangable(hb, d, a, grpd, grpa))
+ {
+ /* Then this hbond/contact will be found again, or it has already been found. */
+ /*return hbNo;*/
+ }
if (bBox){
- if (d>a && bMerge && (bContact || isInterchangable(hb, d, a, grpd, grpa))) { /* acceptor is also a donor and vice versa? */
+ if (d>a && bMerge && isInterchangable(hb, d, a, grpd, grpa)) { /* acceptor is also a donor and vice versa? */
/* return hbNo; */
daSwap = TRUE; /* If so, then their history should be filed with donor and acceptor swapped. */
}
rda2 = iprod(r_da,r_da);
if (bContact) {
- if (daSwap)
+ if (daSwap && grpa == grpd)
return hbNo;
if (rda2 <= rc2){
if (hb->bGem){
if (bDA || (!bDA && (rha2 <= rc2))) {
rvec_sub(x[d],x[hh],r_dh);
if (bBox) {
- if (hb->bGem)
- pbc_correct_gem(r_dh,box,hbox);
- else
- pbc_correct_gem(r_dh,box,hbox);
+ pbc_correct_gem(r_dh,box,hbox);
}
if (!bDA)
fprintf(stderr, "%-7i",data[i]);
}
- static void normalizeACF(real *ct, real *gt, int len)
+ static void normalizeACF(real *ct, real *gt, int nhb, int len)
{
real ct_fac, gt_fac;
int i;
/* Xu and Berne use the same normalization constant */
ct_fac = 1.0/ct[0];
- gt_fac = (gt!=NULL && gt[0]!=0) ? 1.0/gt[0] : 0;
+ gt_fac = (nhb == 0) ? 0 : 1.0/(real)nhb;
+
printf("Normalization for c(t) = %g for gh(t) = %g\n",ct_fac,gt_fac);
for (i=0; i<len; i++)
{
"Ac(t)",
"Cc\\scontact,hb\\v{}\\z{}(t)",
"-dAc\\sfs\\v{}\\z{}/dt" };
- gmx_bool bNorm=FALSE;
+ gmx_bool bNorm=FALSE, bOMP=FALSE;
double nhb = 0;
int nhbi=0;
real *rhbex=NULL,*ht,*gt,*ght,*dght,*kt;
real *ct,*p_ct,tail,tail2,dtail,ct_fac,ght_fac,*cct;
const real tol = 1e-3;
int nframes = hb->nframes,nf;
- unsigned int **h,**g;
+ unsigned int **h=NULL,**g=NULL;
int nh,nhbonds,nhydro,ngh;
t_hbond *hbh;
PSTYPE p, *pfound = NULL, np;
t_E *E;
double *ctdouble, *timedouble, *fittedct;
double fittolerance=0.1;
+ int *dondata=NULL, thisThread;
enum {AC_NONE, AC_NN, AC_GEM, AC_LUZAR};
-
-#ifdef HAVE_OPENMP
- int *dondata=NULL, thisThread;
+#ifdef GMX_OPENMP
+ bOMP = TRUE;
+#else
+ bOMP = FALSE;
#endif
-
printf("Doing autocorrelation ");
/* Decide what kind of ACF calculations to do. */
if (bGemFit)
sprintf(legGem[(bBallistic ? 3:2)], "Ac\\s%s,fit\\v{}\\z{}(t)", gemType);
- } else {
+ }
+ else
+ {
acType = AC_LUZAR;
printf("according to the theory of Luzar and Chandler.\n");
}
nn = nframes/2;
- if (acType != AC_NN ||
-#ifndef HAVE_OPENMP
- TRUE
-#else
- FALSE
-#endif
- ) {
+ if (acType != AC_NN || bOMP) {
snew(h,hb->maxhydro);
snew(g,hb->maxhydro);
}
ngh = 0;
anhb = 0;
- /* ------------------------------------------------
- * I got tired of waiting for the acf calculations
- * and parallelized it with openMP
- * set environment variable CFLAGS = "-fopenmp" when running
- * configure and define DOUSEOPENMP to make use of it.
- */
-
-#ifdef HAVE_OPENMP /* ================================================= \
- * Set up the OpenMP stuff, |
- * like the number of threads and such |
- */
- if (acType != AC_LUZAR)
+ if (acType != AC_LUZAR && bOMP)
{
-/* #if (_OPENMP >= 200805) /\* =====================\ *\/ */
-/* nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
-/* #else */
- nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
-/* #endif /\* _OPENMP >= 200805 ====================/ *\/ */
+ nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
- omp_set_num_threads(nThreads);
+ gmx_omp_set_num_threads(nThreads);
snew(dondata, nThreads);
for (i=0; i<nThreads; i++)
dondata[i] = -1;
fprintf(stderr, "%-7s", tmpstr);
}
}
- fprintf(stderr, "\n"); /* | */
- } /* | */
-#endif /* HAVE_OPENMP ===================================================/ */
+ fprintf(stderr, "\n");
+ }
/* Build the ACF according to acType */
#ifdef HAVE_NN_LOOPS
/* Here we're using the estimated energy for the hydrogen bonds. */
snew(ct,nn);
-#ifdef HAVE_OPENMP /* ==================================\ */
+
#pragma omp parallel \
private(i, j, k, nh, E, rhbex, thisThread), \
default(shared)
{
#pragma omp barrier
- thisThread = omp_get_thread_num();
+ thisThread = gmx_omp_get_thread_num();
rhbex = NULL;
-#endif /* ==============================================/ */
snew(rhbex, n2);
memset(rhbex, 0, n2*sizeof(real)); /* Trust no-one, not even malloc()! */
-#ifdef HAVE_OPENMP /* ################################################## \
- * #
- * #
- */
#pragma omp barrier
#pragma omp for schedule (dynamic)
-#endif
for (i=0; i<hb->d.nrd; i++) /* loop over donors */
{
-#ifdef HAVE_OPENMP /* ====== Write some output ======\ */
+ if (bOMP)
+ {
#pragma omp critical
+ {
+ dondata[thisThread] = i;
+ parallel_print(dondata, nThreads);
+ }
+ }
+ else
{
- dondata[thisThread] = i;
- parallel_print(dondata, nThreads);
+ fprintf(stderr, "\r %i", i);
}
-#else
- fprintf(stderr, "\r %i", i);
-#endif /* ===========================================/ */
for (j=0; j<hb->a.nra; j++) /* loop over acceptors */
{
low_do_autocorr(NULL,oenv,NULL,nframes,1,-1,&(rhbex),hb->time[1]-hb->time[0],
eacNormal,1,FALSE,bNorm,FALSE,0,-1,0,1);
-#ifdef HAVE_OPENMP
#pragma omp critical
-#endif
{
for(k=0; (k<nn); k++)
ct[k] += rhbex[k];
} /* i loop */
sfree(rhbex);
#pragma omp barrier
-#ifdef HAVE_OPENMP
- /* # */
- } /* End of parallel block # */
- /* ##############################################################/ */
- sfree(dondata);
-#endif
+ }
+
+ if (bOMP)
+ {
+ sfree(dondata);
+ }
-
- normalizeACF(ct, NULL, nn);
+ normalizeACF(ct, NULL, 0, nn);
snew(ctdouble, nn);
snew(timedouble, nn);
for (j=0; j<nn; j++)
hb->time[j]-hb->time[0],
ct[j],
ctdouble[j]);
- fclose(fp);
+ xvgrclose(fp);
sfree(ct);
sfree(ctdouble);
sfree(timedouble);
case AC_GEM:
snew(ct,2*n2);
memset(ct,0,2*n2*sizeof(real));
-#ifndef HAVE_OPENMP
+#ifndef GMX_OPENMP
fprintf(stderr, "Donor:\n");
#define __ACDATA ct
#else
#define __ACDATA p_ct
#endif
-#ifdef HAVE_OPENMP /* =========================================\
- * */
-#pragma omp parallel default(none) \
+#pragma omp parallel \
private(i, k, nh, hbh, pHist, h, g, n0, nf, np, j, m, \
pfound, poff, rHbExGem, p, ihb, mMax, \
thisThread, p_ct) \
- shared(hb, dondata, ct, nn, nThreads, n2, stderr, bNorm, \
- nframes, bMerge, bContact)
+ default(shared)
{ /* ########## THE START OF THE ENORMOUS PARALLELIZED BLOCK! ########## */
h = NULL;
g = NULL;
- thisThread = omp_get_thread_num();
+ thisThread = gmx_omp_get_thread_num();
snew(h,hb->maxhydro);
snew(g,hb->maxhydro);
mMax = INT_MIN;
/* I'm using a chunk size of 1, since I expect \
* the overhead to be really small compared \
* to the actual calculations \ */
-#pragma omp for schedule(dynamic,1) nowait /* \ */
-#endif /* HAVE_OPENMP =========================================/ */
-
+#pragma omp for schedule(dynamic,1) nowait
for (i=0; i<hb->d.nrd; i++) {
-#ifdef HAVE_OPENMP
+
+ if (bOMP)
+ {
#pragma omp critical
+ {
+ dondata[thisThread] = i;
+ parallel_print(dondata, nThreads);
+ }
+ }
+ else
{
- dondata[thisThread] = i;
- parallel_print(dondata, nThreads);
+ fprintf(stderr, "\r %i", i);
}
-#else
- fprintf(stderr, "\r %i", i);
-#endif
-
for (k=0; k<hb->a.nra; k++) {
for (nh=0; nh < ((bMerge || bContact) ? 1 : hb->d.nhydro[i]); nh++) {
hbh = hb->hbmap[i][k];
pHist = &(hb->per->pHist[i][k]);
if (ISHB(hbh->history[nh]) && pHist->len != 0) {
-/* No need for a critical section */
-/* #ifdef HAVE_OPENMP */
-/* #pragma omp critical */
-/* #endif */
{
h[nh] = hbh->h[nh];
g[nh] = hb->per->gemtype==gemAD ? hbh->g[nh] : NULL;
srenew(poff,np);
}
-/* This shouldn't have to be critical, right? */
-/* #ifdef HAVE_OPENMP */
-/* #pragma omp critical */
-/* #endif */
{
if (rHbExGem != NULL && rHbExGem[m] != NULL) {
/* This must be done, as this array was most likey
sfree(h);
sfree(g);
-#ifdef HAVE_OPENMP /* =======================================\ */
-#pragma omp critical
+
+ if (bOMP)
{
- for (i=0; i<nn; i++)
- ct[i] += p_ct[i];
+#pragma omp critical
+ {
+ for (i=0; i<nn; i++)
+ ct[i] += p_ct[i];
+ }
+ sfree(p_ct);
}
- sfree(p_ct);
} /* ########## THE END OF THE ENORMOUS PARALLELIZED BLOCK ########## */
- sfree(dondata);
-#endif /* HAVE_OPENMP =======================================/ */
+ if (bOMP)
+ {
+ sfree(dondata);
+ }
- normalizeACF(ct, NULL, nn);
+ normalizeACF(ct, NULL, 0, nn);
fprintf(stderr, "\n\nACF successfully calculated.\n");
fprintf(fp," %10g", fittedct[j]);
fprintf(fp,"\n");
}
- fclose(fp);
+ xvgrclose(fp);
sfree(ctdouble);
sfree(timedouble);
fprintf(stderr,"\n");
sfree(h);
sfree(g);
- normalizeACF(ct, gt, nn);
+ normalizeACF(ct, ght, nhb, nn);
/* Determine tail value for statistics */
tail = 0;
fp = opt2FILE("-hbn",nfile,fnm,"w");
if (opt2bSet("-g",nfile,fnm)) {
fplog = ffopen(opt2fn("-g",nfile,fnm),"w");
- if (bContact)
- fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
- else
- fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
+ fprintf(fplog,"# %10s %12s %12s\n","Donor","Hydrogen","Acceptor");
}
else
fplog = NULL;
ffclose(fplog);
}
-#ifdef HAVE_OPENMP
/* sync_hbdata() updates the parallel t_hbdata p_hb using hb as template.
* It mimics add_frames() and init_frame() to some extent. */
static void sync_hbdata(t_hbdata *hb, t_hbdata *p_hb,
* even though the data its members point to will change,
* hence no need for re-syncing. */
}
-#endif
int gmx_hbond(int argc,char *argv[])
{
"which should contain exactly one atom. In this case, only hydrogen",
"bonds between atoms within the shell distance from the one atom are",
"considered.[PAR]",
+
+ "With option -ac, rate constants for hydrogen bonding can be derived with the model of Luzar and Chandler",
+ "(Nature 394, 1996; J. Chem. Phys. 113:23, 2000) or that of Markovitz and Agmon (J. Chem. Phys 129, 2008).",
+ "If contact kinetics are analyzed by using the -contact option, then",
+ "n(t) can be defined as either all pairs that are not within contact distance r at time t",
+ "(corresponding to leaving the -r2 option at the default value 0) or all pairs that",
+ "are within distance r2 (corresponding to setting a second cut-off value with option -r2).",
+ "See mentioned literature for more details and definitions."
+ "[PAR]",
/* "It is also possible to analyse specific hydrogen bonds with",
"[TT]-sel[tt]. This index file must contain a group of atom triplets",
"Use reversible geminate recombination for the kinetics/thermodynamics calclations. See Markovitch et al., J. Chem. Phys 129, 084505 (2008) for details."},
{ "-diff", FALSE, etREAL, {&D},
"Dffusion coefficient to use in the reversible geminate recombination kinetic model. If negative, then it will be fitted to the ACF along with ka and kd."},
-#ifdef HAVE_OPENMP
+#ifdef GMX_OPENMP
{ "-nthreads", FALSE, etINT, {&nThreads},
"Number of threads used for the parallel loop over autocorrelations. nThreads <= 0 means maximum number of threads. Requires linking with OpenMP. The number of threads is limited by the number of processors (before OpenMP v.3 ) or environment variable OMP_THREAD_LIMIT (OpenMP v.3)"},
#endif
matrix box;
real t,ccut,dist=0.0,ang=0.0;
double max_nhb,aver_nhb,aver_dist;
- int h=0,i,j,k=0,l,start,end,id,ja,ogrp,nsel;
+ int h=0,i=0,j,k=0,l,start,end,id,ja,ogrp,nsel;
int xi,yi,zi,ai;
int xj,yj,zj,aj,xjj,yjj,zjj;
int xk,yk,zk,ak,xkk,ykk,zkk;
gmx_bool bSelected,bHBmap,bStop,bTwo,was,bBox,bTric;
- int *adist,*rdist;
+ int *adist,*rdist,*aptr,*rprt;
int grp,nabin,nrbin,bin,resdist,ihb;
char **leg;
- t_hbdata *hb;
+ t_hbdata *hb,*hbptr;
FILE *fp,*fpins=NULL,*fpnhb=NULL;
t_gridcell ***grid;
t_ncell *icell,*jcell,*kcell;
int threadNr=0;
gmx_bool bGem, bNN, bParallel;
t_gemParams *params=NULL;
+ gmx_bool bEdge_yjj, bEdge_xjj, bOMP;
- CopyRight(stdout,argv[0]);
+ t_hbdata **p_hb=NULL; /* one per thread, then merge after the frame loop */
+ int **p_adist=NULL, **p_rdist=NULL; /* a histogram for each thread. */
+
+#ifdef GMX_OPENMP
+ bOMP = TRUE;
+#else
+ bOMP = FALSE;
+#endif
+
+ CopyRight(stderr,argv[0]);
npargs = asize(pa);
ppa = add_acf_pargs(&npargs,pa);
gmx_fatal(FARGS,"Can not analyze contact between H and A: turn off -noda");
}
}
-
-#ifndef HAVE_LIBGSL
- /* Don't pollute stdout with information about external libraries.
- *
- * printf("NO GSL! Can't find and take away ballistic term in ACF without GSL\n.");
- */
-#endif
/* Initiate main data structure! */
bHBmap = (opt2bSet("-ac",NFILE,fnm) ||
opt2bSet("-hbm",NFILE,fnm) ||
bGem);
-#ifdef HAVE_OPENMP
- /* Same thing here. There is no reason whatsoever to write the specific version of
- * OpenMP used for compilation to stdout for normal usage.
- *
- * printf("Compiled with OpenMP (%i)\n", _OPENMP);
- */
-#endif
-
- /* if (bContact && bGem) */
- /* gmx_fatal(FARGS, "Can't do reversible geminate recombination with -contact yet."); */
-
if (opt2bSet("-nhbdist",NFILE,fnm)) {
const char *leg[MAXHH+1] = { "0 HBs", "1 HB", "2 HBs", "3 HBs", "Total" };
fpnhb = xvgropen(opt2fn("-nhbdist",NFILE,fnm),
bParallel = FALSE;
-#ifndef HAVE_OPENMP
+#ifndef GMX_OPENMP
#define __ADIST adist
#define __RDIST rdist
#define __HBDATA hb
-#else /* HAVE_OPENMP ================================================== \
+#else /* GMX_OPENMP ================================================== \
* Set up the OpenMP stuff, |
* like the number of threads and such |
* Also start the parallel loop. |
#define __ADIST p_adist[threadNr]
#define __RDIST p_rdist[threadNr]
#define __HBDATA p_hb[threadNr]
+#endif
+ if (bOMP)
+ {
+ bParallel = !bSelected;
- bParallel = !bSelected;
+ if (bParallel)
+ {
+ actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, gmx_omp_get_max_threads());
- if (bParallel)
- {
-/* #if (_OPENMP > 200805) */
-/* actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_thread_limit()); */
-/* #else */
- actual_nThreads = min((nThreads <= 0) ? INT_MAX : nThreads, omp_get_num_procs());
-/* #endif */
- omp_set_num_threads(actual_nThreads);
- printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
- fflush(stdout);
- }
- else
- {
- actual_nThreads = 1;
- }
+ gmx_omp_set_num_threads(actual_nThreads);
+ printf("Frame loop parallelized with OpenMP using %i threads.\n", actual_nThreads);
+ fflush(stdout);
+ }
+ else
+ {
+ actual_nThreads = 1;
+ }
- t_hbdata **p_hb; /* one per thread, then merge after the frame loop */
- int **p_adist, **p_rdist; /* a histogram for each thread. */
- snew(p_hb, actual_nThreads);
- snew(p_adist, actual_nThreads);
- snew(p_rdist, actual_nThreads);
- for (i=0; i<actual_nThreads; i++)
- {
- snew(p_hb[i], 1);
- snew(p_adist[i], nabin+1);
- snew(p_rdist[i], nrbin+1);
-
- p_hb[i]->max_frames = 0;
- p_hb[i]->nhb = NULL;
- p_hb[i]->ndist = NULL;
- p_hb[i]->n_bound = NULL;
- p_hb[i]->time = NULL;
- p_hb[i]->nhx = NULL;
-
- p_hb[i]->bHBmap = hb->bHBmap;
- p_hb[i]->bDAnr = hb->bDAnr;
- p_hb[i]->bGem = hb->bGem;
- p_hb[i]->wordlen = hb->wordlen;
- p_hb[i]->nframes = hb->nframes;
- p_hb[i]->maxhydro = hb->maxhydro;
- p_hb[i]->danr = hb->danr;
- p_hb[i]->d = hb->d;
- p_hb[i]->a = hb->a;
- p_hb[i]->hbmap = hb->hbmap;
- p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
- p_hb[i]->per = hb->per;
+ snew(p_hb, actual_nThreads);
+ snew(p_adist, actual_nThreads);
+ snew(p_rdist, actual_nThreads);
+ for (i=0; i<actual_nThreads; i++)
+ {
+ snew(p_hb[i], 1);
+ snew(p_adist[i], nabin+1);
+ snew(p_rdist[i], nrbin+1);
+
+ p_hb[i]->max_frames = 0;
+ p_hb[i]->nhb = NULL;
+ p_hb[i]->ndist = NULL;
+ p_hb[i]->n_bound = NULL;
+ p_hb[i]->time = NULL;
+ p_hb[i]->nhx = NULL;
+
+ p_hb[i]->bHBmap = hb->bHBmap;
+ p_hb[i]->bDAnr = hb->bDAnr;
+ p_hb[i]->bGem = hb->bGem;
+ p_hb[i]->wordlen = hb->wordlen;
+ p_hb[i]->nframes = hb->nframes;
+ p_hb[i]->maxhydro = hb->maxhydro;
+ p_hb[i]->danr = hb->danr;
+ p_hb[i]->d = hb->d;
+ p_hb[i]->a = hb->a;
+ p_hb[i]->hbmap = hb->hbmap;
+ p_hb[i]->time = hb->time; /* This may need re-syncing at every frame. */
+ p_hb[i]->per = hb->per;
#ifdef HAVE_NN_LOOPS
- p_hb[i]->hbE = hb->hbE;
+ p_hb[i]->hbE = hb->hbE;
#endif
- p_hb[i]->nrhb = 0;
- p_hb[i]->nrdist = 0;
+ p_hb[i]->nrhb = 0;
+ p_hb[i]->nrdist = 0;
+ }
}
/* Make a thread pool here,
* instead of forking anew at every frame. */
#pragma omp parallel \
- private(i, j, h, ii, jj, hh, E, \
+ firstprivate(i) \
+ private(j, h, ii, jj, hh, E, \
xi, yi, zi, xj, yj, zj, threadNr, \
dist, ang, peri, icell, jcell, \
grp, ogrp, ai, aj, xjj, yjj, zjj, \
xk, yk, zk, ihb, id, resdist, \
- xkk, ykk, zkk, kcell, ak, k, bTric) \
- default(none) \
- shared(hb, p_hb, p_adist, p_rdist, actual_nThreads, \
- x, bBox, box, hbox, rcut, r2cut, rshell, \
- shatom, ngrid, grid, nframes, t, \
- bParallel, bNN, index, bMerge, bContact, \
- bTwo, bDA,ccut, abin, rbin, top, \
- bSelected, bDebug, stderr, nsel, \
- bGem, oenv, fnm, fpnhb, trrStatus, natoms, \
- status, nabin, nrbin, adist, rdist, debug)
+ xkk, ykk, zkk, kcell, ak, k, bTric, \
+ bEdge_xjj, bEdge_yjj) \
+ default(shared)
{ /* Start of parallel region */
- threadNr = omp_get_thread_num();
-#endif /* HAVE_OPENMP ================================================= */
+ threadNr = gmx_omp_get_thread_num();
+
do
{
+
bTric = bBox && TRICLINIC(box);
-#ifdef HAVE_OPENMP
- sync_hbdata(hb, p_hb[threadNr], nframes, t);
+ if (bOMP)
+ {
+ sync_hbdata(hb, p_hb[threadNr], nframes, t);
+ }
#pragma omp single
-#endif
{
build_grid(hb,x,x[shatom], bBox,box,hbox, (rcut>r2cut)?rcut:r2cut,
rshell, ngrid,grid);
if (hb->bDAnr)
count_da_grid(ngrid, grid, hb->danr[nframes]);
} /* omp single */
-#ifdef HAVE_OPENMP
- p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
-#endif
+
+ if (bOMP)
+ {
+ p_hb[threadNr]->time = hb->time; /* This pointer may have changed. */
+ }
+
if (bNN)
{
#ifdef HAVE_NN_LOOPS /* Unlock this feature when testing */
/* Loop over all atom pairs and estimate interaction energy */
-#ifdef HAVE_OPENMP /* ------- */
+
#pragma omp single
-#endif /* HAVE_OPENMP ------- */
{
addFramesNN(hb, nframes);
}
-#ifdef HAVE_OPENMP /* ---------------- */
+
#pragma omp barrier
#pragma omp for schedule(dynamic)
-#endif /* HAVE_OPENMP ---------------- */
for (i=0; i<hb->d.nrd; i++)
{
for(j=0;j<hb->a.nra; j++)
{
if (bSelected)
{
-#ifdef HAVE_OPENMP
+
#pragma omp single
-#endif
{
/* Do not parallelize this just yet. */
/* int ii; */
} /* if (bSelected) */
else
{
-#ifdef HAVE_OPENMP
+
#pragma omp single
{
-#endif
- if (bGem)
- calcBoxProjection(box, hb->per->P);
+ if (bGem)
+ calcBoxProjection(box, hb->per->P);
+
+ /* loop over all gridcells (xi,yi,zi) */
+ /* Removed confusing macro, DvdS 27/12/98 */
- /* loop over all gridcells (xi,yi,zi) */
- /* Removed confusing macro, DvdS 27/12/98 */
-#ifdef HAVE_OPENMP
}
/* The outer grid loop will have to do for now. */
#pragma omp for schedule(dynamic)
-#endif
for(xi=0; xi<ngrid[XX]; xi++)
for(yi=0; (yi<ngrid[YY]); yi++)
for(zi=0; (zi<ngrid[ZZ]); zi++) {
i = icell->atoms[ai];
/* loop over all adjacent gridcells (xj,yj,zj) */
- /* This is a macro!!! */
- LOOPGRIDINNER(xj,yj,zj,xjj,yjj,zjj,xi,yi,zi,ngrid,bTric) {
- jcell=&(grid[zj][yj][xj].a[ogrp]);
- /* loop over acceptor atoms from other group (ogrp)
- * in this adjacent gridcell (jcell)
- */
- for (aj=0; (aj<jcell->nr); aj++) {
- j = jcell->atoms[aj];
-
- /* check if this once was a h-bond */
- peri = -1;
- ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
- hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
-
- if (ihb) {
- /* add to index if not already there */
- /* Add a hbond */
- add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
-
- /* make angle and distance distributions */
- if (ihb == hbHB && !bContact) {
- if (dist>rcut)
- gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
- ang*=RAD2DEG;
- __ADIST[(int)( ang/abin)]++;
- __RDIST[(int)(dist/rbin)]++;
- if (!bTwo) {
- int id,ia;
- if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
- gmx_fatal(FARGS,"Invalid donor %d",i);
- if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
- gmx_fatal(FARGS,"Invalid acceptor %d",j);
- resdist=abs(top.atoms.atom[i].resind-
- top.atoms.atom[j].resind);
- if (resdist >= max_hx)
- resdist = max_hx-1;
- __HBDATA->nhx[nframes][resdist]++;
+ for(zjj = grid_loop_begin(ngrid[ZZ],zi,bTric,FALSE);
+ zjj <= grid_loop_end(ngrid[ZZ],zi,bTric,FALSE);
+ zjj++)
+ {
+ zj = grid_mod(zjj,ngrid[ZZ]);
+ bEdge_yjj = (zj == 0) || (zj == ngrid[ZZ] - 1);
+ for(yjj = grid_loop_begin(ngrid[YY],yi,bTric,bEdge_yjj);
+ yjj <= grid_loop_end(ngrid[YY],yi,bTric,bEdge_yjj);
+ yjj++)
+ {
+ yj = grid_mod(yjj,ngrid[YY]);
+ bEdge_xjj =
+ (yj == 0) || (yj == ngrid[YY] - 1) ||
+ (zj == 0) || (zj == ngrid[ZZ] - 1);
+ for(xjj = grid_loop_begin(ngrid[XX],xi,bTric,bEdge_xjj);
+ xjj <= grid_loop_end(ngrid[XX],xi,bTric,bEdge_xjj);
+ xjj++)
+ {
+ xj = grid_mod(xjj,ngrid[XX]);
+ jcell=&(grid[zj][yj][xj].a[ogrp]);
+ /* loop over acceptor atoms from other group (ogrp)
+ * in this adjacent gridcell (jcell)
+ */
+ for (aj=0; (aj<jcell->nr); aj++) {
+ j = jcell->atoms[aj];
+
+ /* check if this once was a h-bond */
+ peri = -1;
+ ihb = is_hbond(__HBDATA,grp,ogrp,i,j,rcut,r2cut,ccut,x,bBox,box,
+ hbox,&dist,&ang,bDA,&h,bContact,bMerge,&peri);
+
+ if (ihb) {
+ /* add to index if not already there */
+ /* Add a hbond */
+ add_hbond(__HBDATA,i,j,h,grp,ogrp,nframes,bMerge,ihb,bContact,peri);
+
+ /* make angle and distance distributions */
+ if (ihb == hbHB && !bContact) {
+ if (dist>rcut)
+ gmx_fatal(FARGS,"distance is higher than what is allowed for an hbond: %f",dist);
+ ang*=RAD2DEG;
+ __ADIST[(int)( ang/abin)]++;
+ __RDIST[(int)(dist/rbin)]++;
+ if (!bTwo) {
+ int id,ia;
+ if ((id = donor_index(&hb->d,grp,i)) == NOTSET)
+ gmx_fatal(FARGS,"Invalid donor %d",i);
+ if ((ia = acceptor_index(&hb->a,ogrp,j)) == NOTSET)
+ gmx_fatal(FARGS,"Invalid acceptor %d",j);
+ resdist=abs(top.atoms.atom[i].resind-
+ top.atoms.atom[j].resind);
+ if (resdist >= max_hx)
+ resdist = max_hx-1;
+ __HBDATA->nhx[nframes][resdist]++;
+ }
+ }
+
}
- }
-
- }
- } /* for aj */
- }
- ENDLOOPGRIDINNER;
+ } /* for aj */
+ } /* for xjj */
+ } /* for yjj */
+ } /* for zjj */
} /* for ai */
} /* for grp */
} /* for xi,yi,zi */
} /* if (bSelected) {...} else */
-#ifdef HAVE_OPENMP /* ---------------------------- */
+
/* Better wait for all threads to finnish using x[] before updating it. */
- k = nframes; /* */
-#pragma omp barrier /* */
-#pragma omp critical /* */
- { /* */
+ k = nframes;
+#pragma omp barrier
+#pragma omp critical
+ {
/* Sum up histograms and counts from p_hb[] into hb */
- { /* */
+ if (bOMP)
+ {
hb->nhb[k] += p_hb[threadNr]->nhb[k];
hb->ndist[k] += p_hb[threadNr]->ndist[k];
- for (j=0; j<max_hx; j++) /**/
+ for (j=0; j<max_hx; j++)
hb->nhx[k][j] += p_hb[threadNr]->nhx[k][j];
- } /* */
- } /* */
- /* */
+ }
+ }
+
/* Here are a handful of single constructs
* to share the workload a bit. The most
* important one is of course the last one,
* where there's a potential bottleneck in form
* of slow I/O. */
-#pragma omp single /* ++++++++++++++++, */
-#endif /* HAVE_OPENMP ----------------+------------*/
- { /* + */
- if (hb != NULL) /* */
- { /* + */
+#pragma omp barrier
+#pragma omp single
+ {
+ if (hb != NULL)
+ {
analyse_donor_props(opt2fn_null("-don",NFILE,fnm),hb,k,t,oenv);
- } /* + */
- } /* + */
-#ifdef HAVE_OPENMP /* + */
-#pragma omp single /* +++ +++ */
-#endif /* + */
- { /* + */
- if (fpnhb) /* + */
+ }
+ }
+
+#pragma omp single
+ {
+ if (fpnhb)
do_nhb_dist(fpnhb,hb,t);
- } /* + */
+ }
} /* if (bNN) {...} else + */
-#ifdef HAVE_OPENMP /* + */
-#pragma omp single /* +++ +++ */
-#endif /* + */
- { /* + */
+
+#pragma omp single
+ {
trrStatus = (read_next_x(oenv,status,&t,natoms,x,box));
- nframes++; /* + */
- } /* + */
-#ifdef HAVE_OPENMP /* +++++++++++++++++ */
+ nframes++;
+ }
+
#pragma omp barrier
-#endif
} while (trrStatus);
-#ifdef HAVE_OPENMP
-#pragma omp critical
+ if (bOMP)
{
- hb->nrhb += p_hb[threadNr]->nrhb;
- hb->nrdist += p_hb[threadNr]->nrdist;
- }
- /* Free parallel datastructures */
- sfree(p_hb[threadNr]->nhb);
- sfree(p_hb[threadNr]->ndist);
- sfree(p_hb[threadNr]->nhx);
+#pragma omp critical
+ {
+ hb->nrhb += p_hb[threadNr]->nrhb;
+ hb->nrdist += p_hb[threadNr]->nrdist;
+ }
+ /* Free parallel datastructures */
+ sfree(p_hb[threadNr]->nhb);
+ sfree(p_hb[threadNr]->ndist);
+ sfree(p_hb[threadNr]->nhx);
#pragma omp for
- for (i=0; i<nabin; i++)
- for (j=0; j<actual_nThreads; j++)
+ for (i=0; i<nabin; i++)
+
for (j=0; j<actual_nThreads; j++)
- adist[i] += p_adist[j][i];
+ adist[i] += p_adist[j][i];
#pragma omp for
- for (i=0; i<=nrbin; i++)
- for (j=0; j<actual_nThreads; j++)
- rdist[i] += p_rdist[j][i];
+ for (i=0; i<=nrbin; i++)
+
for (j=0; j<actual_nThreads; j++)
+ rdist[i] += p_rdist[j][i];
- sfree(p_adist[threadNr]);
- sfree(p_rdist[threadNr]);
+ sfree(p_adist[threadNr]);
+ sfree(p_rdist[threadNr]);
+ }
} /* End of parallel region */
- sfree(p_adist);
- sfree(p_rdist);
-#endif
+ if (bOMP)
+ {
+ sfree(p_adist);
+ sfree(p_rdist);
+ }
if(nframes <2 && (opt2bSet("-ac",NFILE,fnm) || opt2bSet("-life",NFILE,fnm)))
{
int id,ia,hh,x,y;
mat.nx=nframes;
- mat.ny=(bContact ? hb->nrdist : hb->nrhb);
+ mat.ny=hb->nrhb;
snew(mat.matrix,mat.nx);
for(x=0; (x<mat.nx); x++)
#include "do_fit.h"
#include "rmpbc.h"
#include "wgms.h"
-#include "magic.h"
#include "pbc.h"
#include "xvgr.h"
#include "xdrf.h"
snew(fnms_out[i],strlen(buf)+32);
sprintf(fnms_out[i],"%d_%s",i,buf);
}
+ sfree(buf);
}
do_demux(nfile_in,fnms,fnms_out,n,val,t,dt_remd,isize,index,dt,oenv);
}
{
searchtime = last_frame_time;
}
- if (xtc_seek_time(stfio,searchtime,fr.natoms))
+ if (xtc_seek_time(stfio,searchtime,fr.natoms,TRUE))
{
gmx_fatal(FARGS,"Error seeking to append position.");
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff