#
# This file is part of the GROMACS molecular simulation package.
#
--# Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
++# Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
# To help us fund GROMACS development, we humbly ask that you cite
# the research papers on the package. Check out http://www.gromacs.org.
-cmake_minimum_required(VERSION 2.8.8)
-# When we require cmake >= 2.8.12, it will provide
-# CMAKE_MINIMUM_REQUIRED_VERSION automatically, but in the meantime we
-# need to set a variable, and it must have a different name.
-set(GMX_CMAKE_MINIMUM_REQUIRED_VERSION "2.8.8")
+cmake_minimum_required(VERSION 3.4.3)
# CMake modules/macros are in a subdirectory to keep this file cleaner
# This needs to be set before project() in order to pick up toolchain files
include(gmxBuildTypeReleaseWithAssert)
if(NOT CMAKE_BUILD_TYPE)
- set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert Profile." FORCE)
- # There's no need to offer a user the choice of ThreadSanitizer
+ set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel Reference RelWithAssert Profile TSAN ASAN MSAN." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release"
- "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert" "Profile")
+ "MinSizeRel" "RelWithDebInfo" "Reference" "RelWithAssert" "Profile" "TSAN" "ASAN" "MSAN")
endif()
if(CMAKE_CONFIGURATION_TYPES)
# Add appropriate GROMACS-specific build types for the Visual
set_property(GLOBAL PROPERTY FIND_LIBRARY_USE_LIB64_PATHS ON)
-# Set a default valgrind suppression file.
-# This unfortunately needs to duplicate information from CTest to work as
-# expected...
-set(MEMORYCHECK_SUPPRESSIONS_FILE
- "${CMAKE_SOURCE_DIR}/cmake/legacy_and_external.supp"
- CACHE FILEPATH
- "File that contains suppressions for the memory checker")
-include(CTest)
+include(gmxCTestUtilities)
+gmx_ctest_init()
include(gmxCPackUtilities)
gmx_cpack_init()
option(GMX_USE_OPENCL "Enable OpenCL acceleration" OFF)
-# Decide on GPU settings based on user-settings and GPU/CUDA detection.
-# GCC 4.6 requires CUDA 5.0 and VS2015 requires CUDA 8.0
+# Decide on GPU settings based on user-settings and GPU/CUDA
+# detection. GCC 4.8 requires CUDA 6.0 (but we choose 6.5 for the
+# preliminary C++11 support), icc 15 requires CUDA 7.0, and VS2015
+# requires CUDA 8.0
if(MSVC)
set(REQUIRED_CUDA_VERSION 8.0)
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "Intel")
+ set(REQUIRED_CUDA_VERSION 7.0)
else()
- set(REQUIRED_CUDA_VERSION 5.0)
+ set(REQUIRED_CUDA_VERSION 6.5)
endif()
set(REQUIRED_CUDA_COMPUTE_CAPABILITY 2.0)
option(GMX_BROKEN_CALLOC "Work around broken calloc()" OFF)
mark_as_advanced(GMX_BROKEN_CALLOC)
- option(GMX_LOAD_PLUGINS "Compile with plugin support, needed to read VMD supported file formats" ON)
- mark_as_advanced(GMX_LOAD_PLUGINS)
option(GMX_OPENMP "Enable OpenMP-based multithreading" ON)
include(gmxCFlags)
gmx_c_flags()
+# These collect libraries that Gromacs requires for linking:
+# - Libraries that are required for libgromacs (only)
+set(GMX_EXTRA_LIBRARIES "")
+# - Libraries that are required for all code in the repository
+set(GMX_COMMON_LIBRARIES "")
+# - Libraries that all code linked against libgromacs needs
+# (i.e., something that is exposed in installed headers).
+set(GMX_PUBLIC_LIBRARIES "")
+
# This variable should be used for additional compiler flags which are not
# generated in gmxCFlags nor are SIMD or MPI related.
set(EXTRA_C_FLAGS "")
# Run through a number of tests for buggy compilers and other issues
include(gmxTestCompilerProblems)
gmx_test_compiler_problems()
-# GMX_SIMD will not be set automatically until the second
-# pass (which is not strictly guaranteed to occur), so putting this
-# check here among logically-related tests is inefficient, but the
-# potential loss is likely zero.
-if(GMX_SIMD STREQUAL "AVX_256"
- AND CMAKE_COMPILER_IS_GNUCC
- AND (C_COMPILER_VERSION VERSION_EQUAL "4.6.1"
- OR CXX_COMPILER_VERSION VERSION_EQUAL "4.6.1"))
- message(FATAL_ERROR "gcc 4.6.1 has buggy support for AVX, and GROMACS mdrun will not work. If you want simulation performance, use a more recent compiler. Otherwise, use GMX_SIMD=SSE4.1")
- # See http://gcc.gnu.org/bugzilla/show_bug.cgi?id=49002
-endif()
# Implement double-precision option. This is complicated because we
# need installed headers to use the precision mode of the build that
# Find external packages #
########################################################################
-# TNG wants zlib if it is available. And static libxml2 might have a dependency
-find_package(ZLIB QUIET)
-include(gmxTestZLib)
-gmx_test_zlib(HAVE_ZLIB)
-
# Unconditionally find the package, as it is also required for unit
# tests. This exports LIBXML2_FOUND, which we should not use because
# it does not tell us that linking will succeed. Instead, we test that
option(GMX_HWLOC "Add support for hwloc Portable Hardware locality library" ${GMX_HWLOC_DEFAULT})
if(GMX_HWLOC)
if(HWLOC_FOUND)
- include_directories(${HWLOC_INCLUDE_DIRS})
+ include_directories(SYSTEM ${HWLOC_INCLUDE_DIRS})
list(APPEND GMX_EXTRA_LIBRARIES ${HWLOC_LIBRARIES})
else()
message(FATAL_ERROR "Hwloc package support requested, but not found.")
add_definitions(-DTMPI_ATOMICS_DISABLED)
endif()
-# Note this relies on zlib detection having already run
include(gmxManageTNG)
include(gmxManageLmfit)
if(GMX_GPU)
# now that we have detected the dependencies, do the second configure pass
gmx_gpu_setup()
-else()
- mark_as_advanced(CUDA_HOST_COMPILER)
endif()
if(CYGWIN)
add_definitions( -DHAVE_CONFIG_H )
include_directories(BEFORE ${CMAKE_SOURCE_DIR}/src)
-include_directories(BEFORE ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/include)
+# TODO required at high level because both libgromacs and progs/mdrun
+# require it, both for thread-MPI and its atomics and mutexes.
+include_directories(BEFORE SYSTEM ${CMAKE_SOURCE_DIR}/src/external/thread_mpi/include)
# Required for config.h, maybe should only be set in src/CMakeLists.txt
include_directories(BEFORE ${CMAKE_BINARY_DIR}/src)
include(gmxManageLinearAlgebraLibraries)
- # Whether GROMACS will really try to compile support for VMD plugins
- set(GMX_USE_PLUGINS OFF)
+ include(gmxManagePluginSupport)
- if(GMX_LOAD_PLUGINS)
- if(NOT WIN32)
- # Native Windows does not have, nor need dlopen
- include(gmxTestdlopen)
- gmx_test_dlopen(HAVE_DLOPEN)
- endif()
-
- # so, should we use plug-ins?
- if(WIN32 OR (HAVE_DLOPEN AND BUILD_SHARED_LIBS))
- if(NOT VMD_QUIETLY)
- MESSAGE(STATUS "Using dynamic plugins (e.g VMD-supported file formats)")
- endif()
+ if (GMX_USE_PLUGINS)
if(NOT GMX_VMD_PLUGIN_PATH)
- find_package(VMD)
+ find_package(VMD)
endif()
- 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()
- set(VMD_QUIETLY TRUE CACHE INTERNAL "")
# Link real-time library for POSIX timers. The check for clock_gettime
# confirms the linkability of rt.
endif()
# Math and thread libraries must often come after all others when linking...
-if(HAVE_LIBM)
- list(APPEND GMX_EXTRA_LIBRARIES m)
+if (HAVE_LIBM)
+ list(APPEND GMX_PUBLIC_LIBRARIES m)
endif()
option(GMX_NACL "Configure for Native Client builds" OFF)
message("CMAKE_CXX_FLAGS_RELEASE: ${GMXC_CXXFLAGS_RELEASE}")
message("CMAKE_CXX_FLAGS_DEBUG: ${GMXC_CXXFLAGS_DEBUG}")
message("CMAKE_EXE_LINKER_FLAGS: ${FFT_LINKER_FLAGS} ${MPI_LINKER_FLAGS}")
- message("CMAKE_SHARED_LINKER_FLAGS: ${MPI_LINKER_FLAGS}")
+ message("CMAKE_SHARED_LINKER_FLAGS: ${FFT_LINKER_FLAGS} ${MPI_LINKER_FLAGS}")
endif()
if(NOT GMX_OPENMP)
################################################################
# Shared library load path settings
################################################################
-# CMake supports RPATH on OS X only from 2.8.12 upwards.
-# CMAKE_SYSTEM_VERSION > 8.0 matches OS X 10.5 and above, where RPATH support
-# was added.
-
if(NOT GMX_BUILD_SHARED_EXE)
# No rpath
set(CMAKE_SKIP_RPATH TRUE)
set(CMAKE_EXE_LINK_DYNAMIC_C_FLAGS) # remove -Wl,-Bdynamic
set(CMAKE_EXE_LINK_DYNAMIC_CXX_FLAGS)
-elseif((NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin") OR
- ((CMAKE_SYSTEM_VERSION VERSION_GREATER 8.0) AND (NOT CMAKE_VERSION VERSION_LESS 2.8.12)))
+else()
# The build folder always has bin/ and lib/; if we are also going to
# install to lib/, then the installation RPATH works also in the build
# tree. This makes installation slightly faster (no need to rewrite the
endif()
# Set the RPATH as relative to the executable location to make the
# binaries relocatable.
- if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin")
+ if(NOT CMAKE_SYSTEM_NAME STREQUAL "Darwin") #Assume OS X >=10.5
set(CMAKE_INSTALL_RPATH "\$ORIGIN/../${GMX_LIB_INSTALL_DIR}")
else()
set(CMAKE_INSTALL_RPATH "@executable_path/../${GMX_LIB_INSTALL_DIR}")
endif()
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
set(CMAKE_MACOSX_RPATH 1)
-else()
- # We are on Darwin/OSX, and CMake cannot handle RPATHs automatically.
- if(CMAKE_SYSTEM_VERSION VERSION_GREATER 8.0)
- # Set the RPATH options manually.
- set(CMAKE_INSTALL_NAME_DIR "@rpath")
- set(GMX_EXE_LINKER_FLAGS ${GMX_EXE_LINKER_FLAGS} "-Wl,-rpath,@executable_path/../${GMX_LIB_INSTALL_DIR}")
- else()
- # Use the old INSTALL_NAME_DIR mechanism if RPATH is not supported.
- set(CMAKE_INSTALL_NAME_DIR "${CMAKE_INSTALL_PREFIX}/${LIB_INSTALL_DIR}")
- endif()
endif()
#COPYING file: Only necessary for binary distributions.
endif()
if (BUILD_TESTING)
- # "tests" target builds all the separate test binaries.
- add_custom_target(tests)
- # "run-ctest" is an internal target that actually runs the tests.
- # This is necessary to be able to add separate targets that execute as part
- # of 'make check', but are ensured to be executed after the actual tests.
- add_custom_target(run-ctest
- COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
- COMMENT "Running all tests"
- VERBATIM)
- add_dependencies(run-ctest tests)
- # "check" target builds and runs all tests.
- add_custom_target(check DEPENDS run-ctest)
+ include(tests/CheckTarget.cmake)
endif()
if (NOT GMX_BUILD_MDRUN_ONLY)
#
# This file is part of the GROMACS molecular simulation package.
#
- # Copyright (c) 2012,2013,2014,2015, by the GROMACS development team, led by
-# Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
++# Copyright (c) 2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
endif()
#Verify FFTW is compiled with fPIC (necessary for shared libraries)
- if (CMAKE_OBJDUMP AND CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" AND BUILD_SHARED_LIBS AND NOT CYGWIN)
+ if (CMAKE_OBJDUMP AND CMAKE_SYSTEM_PROCESSOR STREQUAL "x86_64" AND BUILD_SHARED_LIBS AND NOT CYGWIN AND NOT APPLE)
execute_process(COMMAND ${CMAKE_OBJDUMP} --reloc ${${FFTW}_LIBRARY} OUTPUT_VARIABLE ${FFTW}_OBJDUMP)
if (${${FFTW}_OBJDUMP} MATCHES "R_X86_64" #Should always be true for static libraries. Checks that objdump works properly and that the library isn't dynamic
AND NOT ${${FFTW}_OBJDUMP} MATCHES "R_X86_64_PLT32")
#
# This file is part of the GROMACS molecular simulation package.
#
- # Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
+ # Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
# Problematic with CUDA
# GMX_TEST_CXXFLAG(CXXFLAGS_WARN_EFFCXX "-Wnon-virtual-dtor" GMXC_CXXFLAGS)
GMX_TEST_CXXFLAG(CXXFLAGS_WARN_EXTRA "-Wextra -Wno-missing-field-initializers -Wpointer-arith" GMXC_CXXFLAGS)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN_UNDEF "-Wundef" GMXC_CXXFLAGS)
+ # CUDA versions prior to 7.5 come with a header (math_functions.h) which uses the _MSC_VER macro
+ # unconditionally, so we don't use -Wundef for earlier CUDA versions.
+ if(NOT(GMX_GPU AND CUDA_VERSION VERSION_LESS "7.5"))
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN_UNDEF "-Wundef" GMXC_CXXFLAGS)
+ endif()
GMX_TEST_CFLAG(CXXFLAGS_WARN_REL "-Wno-array-bounds" GMXC_CXXFLAGS_RELEASE_ONLY)
endif()
# new in gcc 4.5
if (CMAKE_C_COMPILER_ID MATCHES "Intel")
if (NOT WIN32)
if(NOT GMX_OPENMP)
- if(CMAKE_C_COMPILER_VERSION VERSION_GREATER 13.99.99)
# 3180: unrecognized OpenMP #pragma
- GMX_TEST_CFLAG(CFLAGS_PRAGMA "-wd3180" GMXC_CFLAGS)
- else()
-# 161: unrecognized #pragma
- GMX_TEST_CFLAG(CFLAGS_PRAGMA "-wd161" GMXC_CFLAGS)
- endif()
+ GMX_TEST_CFLAG(CFLAGS_PRAGMA "-wd3180" GMXC_CFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- if(CMAKE_C_COMPILER_VERSION VERSION_LESS 15.00.00)
-# 193: zero used for undefined preprocessing identifier ".."
- GMX_TEST_CFLAG(CFLAGS_WARN_OLD -wd193 GMXC_CFLAGS)
- endif()
# 177: function/variable ".." was declared but never referenced
-# 271: trailing comma is nonstandard
-# 304: access control not specified ("public" by default)
-# 383: value copied to temporary, reference to temporary used
-# 424: extra ";" ignored
-# 444: destructor for base class ".." is not virtual
-# 522: function ".." redeclared "inline" after being called
+# 411: class defines no constructor to initialize the following (incorrect for struct, initializer list works)
# 593: variable ".." was set but never used
-# 869: parameter ".." was never referenced
# 981: operands are evaluated in unspecified order
#1418: external function definition with no prior declaration
#1419: external declaration in primary source file
#2547: ".." was specified as both a system and non-system include directory
#2557: comparison between signed and unsigned operands
#3280: declaration hides member ".."
-#3346: dynamic exception specifications are deprecated
#11074: Inlining inhibited by limit max-size(/max-total-size)
#11076: To get full report use -opt-report=3 -opt-report-phase ipo (shown for previous remark)
- GMX_TEST_CFLAG(CFLAGS_WARN "-w3 -wd177 -wd271 -wd304 -wd383 -wd424 -wd444 -wd522 -wd593 -wd869 -wd981 -wd1418 -wd1419 -wd1572 -wd1599 -wd2259 -wd2415 -wd2547 -wd2557 -wd3280 -wd3346 -wd11074 -wd11076" GMXC_CFLAGS)
+ GMX_TEST_CFLAG(CFLAGS_WARN "-w3 -wd177 -wd411 -wd593 -wd981 -wd1418 -wd1419 -wd1572 -wd1599 -wd2259 -wd2415 -wd2547 -wd2557 -wd3280 -wd11074 -wd11076" GMXC_CFLAGS)
endif()
GMX_TEST_CFLAG(CFLAGS_STDGNU "-std=gnu99" GMXC_CFLAGS)
- GMX_TEST_CFLAG(CFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias" GMXC_CFLAGS_RELEASE)
+ GMX_TEST_CFLAG(CFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -static-intel" GMXC_CFLAGS_RELEASE)
GMX_TEST_CFLAG(CFLAGS_DEBUG "-O0" GMXC_CFLAGS_DEBUG) #icc defaults to -O2 even with -g
GMX_TEST_CFLAG(CFLAGS_FP_RELASSERT "-fp-model except -fp-model precise" GMXC_CFLAGS_RELWITHASSERT)
else()
if(NOT GMX_OPENMP)
- if(CMAKE_C_COMPILER_VERSION VERSION_GREATER 13.99.99)
- GMX_TEST_CFLAG(CFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
- else()
- GMX_TEST_CFLAG(CFLAGS_PRAGMA "/wd161" GMXC_CFLAGS)
- endif()
+ GMX_TEST_CFLAG(CFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- if(CMAKE_C_COMPILER_VERSION VERSION_LESS 15.00.00)
- GMX_TEST_CFLAG(CFLAGS_WARN_OLD /wd193 GMXC_CFLAGS)
- endif()
- GMX_TEST_CFLAG(CFLAGS_WARN "/W3 /wd177 /wd271 /wd304 /wd383 /wd424 /wd444 /wd522 /wd593 /wd869 /wd981 /wd1418 /wd1419 /wd1572 /wd1599 /wd2259 /wd2415 /wd2547 /wd2557 /wd3280 /wd3346" GMXC_CFLAGS)
+GMX_TEST_CFLAG(CFLAGS_WARN "/W3 /wd177 /wd411 /wd593 /wd981 /wd1418 /wd1419 /wd1572 /wd1599 /wd2259 /wd2415 /wd2547 /wd2557 /wd3280" GMXC_CFLAGS)
endif()
GMX_TEST_CFLAG(CFLAGS_OPT "/Qip" GMXC_CFLAGS_RELEASE)
endif()
if (CMAKE_CXX_COMPILER_ID MATCHES "Intel")
if (NOT WIN32)
if(NOT GMX_OPENMP)
- if(CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 13.99.99)
- GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "-wd3180" GMXC_CXXFLAGS)
- else()
- GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "-wd161" GMXC_CXXFLAGS)
- endif()
+ GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "-wd3180" GMXC_CXXFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 15.00.00)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN_OLD -wd193 GMXC_CXXFLAGS)
+ if (GMX_GPU)
+# Suppress warnings from CUDA headers
+# 7: unrecognized token
+# 82: storage class is not first
+# The below are also required for math_functions.h / math_functions.hpp at least until CUDA 8.0-RC
+# 193: zero used for undefined preprocessing identifer
+# 3346:dynamic exception specifiers are deprecated
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN_OLD_GPU "-wd7 -wd82 -wd193 -wd3346" GMXC_CXXFLAGS)
endif()
#All but the following warnings are identical for the C-compiler (see above)
-#1782: #pragma once is obsolete
+# 383: value copied to temporary, reference to temporary used
+# 444: destructor for base class ".." is not virtual
#2282: unrecognized GCC pragma
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "-w3 -wd177 -wd271 -wd304 -wd383 -wd424 -wd444 -wd522 -wd593 -wd869 -wd981 -wd1418 -wd1419 -wd1572 -wd1599 -wd2259 -wd2415 -wd2547 -wd2557 -wd3280 -wd3346 -wd11074 -wd11076 -wd1782 -wd2282" GMXC_CXXFLAGS)
+ GMX_TEST_CXXFLAG(CXXFLAGS_WARN "-w3 -wd177 -wd383 -wd411 -wd444 -wd981 -wd1418 -wd1572 -wd1599 -wd2259 -wd3280 -wd11074 -wd11076 -wd2282" GMXC_CXXFLAGS)
endif()
- GMX_TEST_CXXFLAG(CXXFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias" GMXC_CXXFLAGS_RELEASE)
+ GMX_TEST_CXXFLAG(CXXFLAGS_OPT "-ip -funroll-all-loops -alias-const -ansi-alias -no-prec-div -fimf-domain-exclusion=14 -static-intel" GMXC_CXXFLAGS_RELEASE)
GMX_TEST_CXXFLAG(CXXFLAGS_DEBUG "-O0" GMXC_CXXFLAGS_DEBUG)
GMX_TEST_CXXFLAG(CXXFLAGS_FP_RELASSERT "-fp-model except -fp-model precise" GMXC_CXXFLAGS_RELWITHASSERT)
else()
- #809: exception specification for virtual function X is incompatible with that of overridden function
if(NOT GMX_OPENMP)
- if(CMAKE_CXX_COMPILER_VERSION VERSION_GREATER 13.99.99)
- GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
- else()
- GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "/wd161" GMXC_CXXFLAGS)
- endif()
+ GMX_TEST_CXXFLAG(CXXFLAGS_PRAGMA "/wd3180" GMXC_CFLAGS)
endif()
if (GMX_COMPILER_WARNINGS)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/W3 /wd177 /wd383 /wd411 /wd444 /wd809 /wd981 /wd1418 /wd1572 /wd1599 /wd2259 /wd3280 /wd11074 /wd11076 /wd2282" GMXC_CXXFLAGS)
- if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS 15.00.00)
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN_OLD /wd193 GMXC_CXXFLAGS)
- endif()
+ #809: exception specification for virtual function X is incompatible with that of overridden function
- GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/W3 /wd177 /wd271 /wd304 /wd383 /wd424 /wd444 /wd522 /wd593 /wd809 /wd869 /wd981 /wd1418 /wd1419 /wd1572 /wd1599 /wd2259 /wd2415 /wd2547 /wd2557 /wd3280 /wd3346 /wd1782 /wd2282" GMXC_CXXFLAGS)
++ GMX_TEST_CXXFLAG(CXXFLAGS_WARN "/W3 /wd177 /wd383 /wd411 /wd444 /wd809 /wd981 /wd1418 /wd1572 /wd1599 /wd2259 /wd3280 /wd11074 /wd11076 /wd2282" GMXC_CXXFLAGS)
endif()
GMX_TEST_CXXFLAG(CXXFLAGS_OPT "/Qip" GMXC_CXXFLAGS_RELEASE)
endif()
--- /dev/null
-# Copyright (c) 2016, by the GROMACS development team, led by
+ #
+ # This file is part of the GROMACS molecular simulation package.
+ #
++# Copyright (c) 2016,2017, by the GROMACS development team, led by
+ # Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ # and including many others, as listed in the AUTHORS file in the
+ # top-level source directory and at http://www.gromacs.org.
+ #
+ # GROMACS is free software; you can redistribute it and/or
+ # modify it under the terms of the GNU Lesser General Public License
+ # as published by the Free Software Foundation; either version 2.1
+ # of the License, or (at your option) any later version.
+ #
+ # GROMACS is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ # Lesser General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public
+ # License along with GROMACS; if not, see
+ # http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ #
+ # If you want to redistribute modifications to GROMACS, please
+ # consider that scientific software is very special. Version
+ # control is crucial - bugs must be traceable. We will be happy to
+ # consider code for inclusion in the official distribution, but
+ # derived work must not be called official GROMACS. Details are found
+ # in the README & COPYING files - if they are missing, get the
+ # official version at http://www.gromacs.org.
+ #
+ # To help us fund GROMACS development, we humbly ask that you cite
+ # the research papers on the package. Check out http://www.gromacs.org.
+
+ include(gmxOptionUtilities)
+
+ # Sets GMX_USE_PLUGINS=ON in the parent scope if the toolchain and
+ # user selections permit the build to support plugin loading.
+ function(gmx_manage_plugin_support)
+ gmx_option_trivalue(GMX_LOAD_PLUGINS "Compile with plugin support, needed to read VMD supported file formats" AUTO)
+ mark_as_advanced(GMX_LOAD_PLUGINS)
+
+ # Find out if non-Windows builds can support plugins. Native Windows
+ # neither needs nor has library support.
+ if (NOT WIN32)
+ # TODO Make a proper find_package for dlopen to find
+ # dlfcn.h. The CMake variable CMAKE_DL_LIBS works magically
+ # for the library, however.
+ include(gmxTestdlopen)
+ gmx_test_dlopen(HAVE_DLOPEN)
+ endif()
+
+ # Keep the status line quiet unless something relevant has
+ # changed.
+ gmx_check_if_changed(EMIT_STATUS_MESSAGES GMX_LOAD_PLUGINS BUILD_SHARED_LIBS HAVE_DLOPEN)
+
+ # Whether GROMACS will really try to compile support for VMD
+ # plugins.
+ set(GMX_USE_PLUGINS OFF)
+
+ # Plugins are supported natively on Windows
+ if (WIN32 OR (BUILD_SHARED_LIBS AND HAVE_DLOPEN))
+ set(GMX_USE_PLUGINS ${GMX_LOAD_PLUGINS})
+ elseif(GMX_LOAD_PLUGINS)
+ # Can't support plugins for some reason. If the user required
+ # plugins, emit fatal errors. Otherwise, emit status messages
+ # for AUTO and be silent for OFF.
+ set(message "")
+ if (NOT HAVE_DLOPEN)
+ set(message "${message}dlopen() support for using dynamic plugins for VMD-supported file formats is missing. ")
+ endif()
+ if(NOT BUILD_SHARED_LIBS)
+ set(message "${message}GROMACS only supports plugins in a build that uses shared libraries, which can be disabled for various reasons. BUILD_SHARED_LIBS=on and a toolchain that supports dynamic linking is required. (Hint: GMX_PREFER_STATIC_LIBS and GMX_BUILD_MDRUN_ONLY can influence the default BUILD_SHARED_LIBS, so if you need plugins, reconsider those choices.) ")
+ endif()
+ if (GMX_LOAD_PLUGINS_FORCE)
+ message(FATAL_ERROR "${message}Cannot build with GMX_LOAD_PLUGINS=${GMX_LOAD_PLUGINS}.")
+ endif()
+ if (GMX_LOAD_PLUGINS_AUTO AND EMIT_STATUS_MESSAGES)
+ message(STATUS "${message}")
+ endif()
+ endif()
+
+ if(EMIT_STATUS_MESSAGES)
+ if(GMX_USE_PLUGINS)
+ MESSAGE(STATUS "Using dynamic plugins (e.g VMD-supported file formats)")
+ else()
+ MESSAGE(STATUS "Not using dynamic plugins (e.g VMD-supported file formats)")
+ endif()
+ endif()
+ set(GMX_USE_PLUGINS ${GMX_USE_PLUGINS} PARENT_SCOPE)
+ endfunction()
+
+ gmx_manage_plugin_support()
+
+ if(GMX_USE_PLUGINS)
+ 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}")
+ endif()
#
# This file is part of the GROMACS molecular simulation package.
#
- # Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by
+ # Copyright (c) 2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
# First check that the compiler is OK, and find the appropriate flag.
if(WIN32 AND NOT MINGW)
- set(CXX11_CXX_FLAG "/Qstd=c++0x")
+ set(CXX11_CXX_FLAG "/Qstd=c++11")
elseif(CYGWIN)
- set(CXX11_CXX_FLAG "-std=gnu++0x") #required for strdup
+ set(CXX11_CXX_FLAG "-std=gnu++11") #required for strdup
else()
- set(CXX11_CXX_FLAG "-std=c++0x")
+ set(CXX11_CXX_FLAG "-std=c++11")
endif()
CHECK_CXX_COMPILER_FLAG("${CXX11_CXX_FLAG}" CXXFLAG_STD_CXX0X)
if(NOT CXXFLAG_STD_CXX0X)
int main() {
// Test nullptr
double *x = nullptr;
+ (void)x; // Suppressing unused variable warning
// Test range-based for loops
int array[5] = { 1, 2, 3, 4, 5 };
for (int& x : array)
x *= 2;
+ // Test alignas
+ alignas(4*sizeof(int)) int y;
}" CXX11_SUPPORTED)
+ if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+ if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.8.1")
+ message(FATAL_ERROR "GROMACS requires version 4.8.1 or later of the GNU C++ compiler for complete C++11 support")
+ endif()
+ elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
+ if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "3.3")
+ message(FATAL_ERROR "GROMACS requires version 3.3 or later of the Clang C++ compiler for complete C++11 support")
+ endif()
+ elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+ if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "15.0")
+ message(FATAL_ERROR "GROMACS requires version 15.0 or later of the Intel C++ compiler for complete C++11 support")
+ endif()
+ elseif("${CMAKE_CXX_COMPILER_ID}" STREQUAL "MSVC")
+ if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "19.0.23026")
+ message(FATAL_ERROR "GROMACS requires version 2015 (19.0.23026) or later of the MSVC C++ compiler for complete C++11 support")
+ endif()
+ endif()
if(CXX11_SUPPORTED)
set(${CXX11_CXX_FLAG_NAME} ${CXX11_CXX_FLAG} PARENT_SCOPE)
else()
intPointer p(new int(10));
std::map<int, std::unique_ptr<int>> m;
m.insert(std::make_pair(5, std::move(p)));
- auto start = std::chrono::system_clock::now();
- if (std::chrono::system_clock::now() - start < std::chrono::seconds(2))
+ auto start = std::chrono::steady_clock::now();
+ if (std::chrono::steady_clock::now() - start < std::chrono::seconds(2))
{
std::thread t;
}
#include <stdlib.h>
#include <string.h>
+#include <cassert>
+
#include <algorithm>
#include <string>
print_missing_interactions_atoms(fplog, cr, top_global, &top_local->idef);
write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
- -1, state_local->x, state_local->box);
+ -1, as_rvec_array(state_local->x.data()), state_local->box);
std::string errorMessage;
}
}
-/*! \brief Run the reverse ilist generation and store it when \p bAssign = TRUE */
+/*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */
static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom,
const int * const * vsite_pbc,
int *count,
a = ia[1+link];
if (bAssign)
{
+ assert(r_il); //with bAssign not allowed to be null
+ assert(r_index);
r_il[r_index[a]+count[a]] =
(ftype == F_CONSTRNC ? F_CONSTR : ftype);
r_il[r_index[a]+count[a]+1] = ia[0];
snew(count, nat_mt);
low_make_reverse_ilist(ilist, atoms->atom, vsite_pbc,
count,
- bConstr, bSettle, bBCheck, NULL, NULL,
+ bConstr, bSettle, bBCheck, nullptr, nullptr,
bLinkToAllAtoms, FALSE);
snew(ril_mt->index, nat_mt+1);
/* Make the atom to interaction list for this molecule type */
nint_mt[mt] =
make_reverse_ilist(molt->ilist, &molt->atoms,
- vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
+ vsite_pbc_molt ? vsite_pbc_molt[mt] : nullptr,
rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
&rt->ril_mt[mt]);
{
t_atoms atoms_global;
- rt->ril_intermol.index = NULL;
- rt->ril_intermol.il = NULL;
+ rt->ril_intermol.index = nullptr;
+ rt->ril_intermol.il = nullptr;
atoms_global.nr = mtop->natoms;
- atoms_global.atom = NULL; /* Only used with virtual sites */
+ atoms_global.atom = nullptr; /* Only used with virtual sites */
*nint +=
make_reverse_ilist(mtop->intermolecular_ilist, &atoms_global,
- NULL,
+ nullptr,
rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
&rt->ril_intermol);
}
rt->nthread = gmx_omp_nthreads_get(emntDomdec);
snew(rt->th_work, rt->nthread);
- if (vsite_pbc_molt != NULL)
+ if (vsite_pbc_molt != nullptr)
{
for (thread = 0; thread < rt->nthread; thread++)
{
*/
dd->reverse_top = make_reverse_top(mtop, ir->efep != efepNO,
- vsite ? vsite->vsite_pbc_molt : NULL,
+ vsite ? vsite->vsite_pbc_molt : nullptr,
!dd->bInterCGcons, !dd->bInterCGsettles,
bBCheck, &dd->nbonded_global);
int nral1 = 0, ftv = 0;
vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
- vsite->vsite_pbc_loc != NULL);
+ vsite->vsite_pbc_loc != nullptr);
if (vpbc)
{
nral1 = 1 + NRAL(ftype);
/* Position restraints need an additional treatment */
if (ftype == F_POSRES || ftype == F_FBPOSRES)
{
- int nposres = dest->il[ftype].nr/2;
- if (nposres > dest->iparams_posres_nalloc)
+ int nposres = dest->il[ftype].nr/2;
+ // TODO: Simplify this code using std::vector
+ t_iparams * &iparams_dest = (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
+ int &posres_nalloc = (ftype == F_POSRES ? dest->iparams_posres_nalloc : dest->iparams_fbposres_nalloc);
+ if (nposres > posres_nalloc)
{
- dest->iparams_posres_nalloc = over_alloc_large(nposres);
- srenew(dest->iparams_posres, dest->iparams_posres_nalloc);
+ posres_nalloc = over_alloc_large(nposres);
+ srenew(iparams_dest, posres_nalloc);
}
+
/* Set nposres to the number of original position restraints in dest */
for (int s = 1; s < nsrc; s++)
{
nposres -= src[s].idef.il[ftype].nr/2;
}
+
for (int s = 1; s < nsrc; s++)
{
- if (ftype == F_POSRES)
- {
- for (int i = 0; i < src[s].idef.il[ftype].nr/2; i++)
- {
- /* Correct the index into iparams_posres */
- dest->il[ftype].iatoms[nposres*2] = nposres;
- /* Copy the position restraint force parameters */
- dest->iparams_posres[nposres] = src[s].idef.iparams_posres[i];
- nposres++;
- }
- }
- else
+ const t_iparams *iparams_src = (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
+
+ for (int i = 0; i < src[s].idef.il[ftype].nr/2; i++)
{
- for (int i = 0; i < src[s].idef.il[ftype].nr/2; i++)
- {
- /* Correct the index into iparams_fbposres */
- dest->il[ftype].iatoms[nposres*2] = nposres;
- /* Copy the position restraint force parameters */
- dest->iparams_fbposres[nposres] = src[s].idef.iparams_fbposres[i];
- nposres++;
- }
+ /* Correct the index into iparams_posres */
+ dest->il[ftype].iatoms[nposres*2] = nposres;
+ /* Copy the position restraint force parameters */
+ iparams_dest[nposres] = iparams_src[i];
+ nposres++;
}
}
}
}
else
{
- vsite_pbc = NULL;
- vsite_pbc_nalloc = NULL;
+ vsite_pbc = nullptr;
+ vsite_pbc_nalloc = nullptr;
}
rt->th_work[thread].nbonded =
real rc = -1;
ivec rcheck;
int d, nexcl;
- t_pbc pbc, *pbc_null = NULL;
+ t_pbc pbc, *pbc_null = nullptr;
if (debug)
{
}
else
{
- pbc_null = NULL;
+ pbc_null = nullptr;
}
}
}
for (i = 0; i < F_NRE; i++)
{
- top->idef.il[i].iatoms = NULL;
+ top->idef.il[i].iatoms = nullptr;
top->idef.il[i].nalloc = 0;
}
top->idef.ilsort = ilsortUNKNOWN;
buf[1] = state_global->ngtc;
buf[2] = state_global->nnhpres;
buf[3] = state_global->nhchainlength;
- buf[4] = state_global->dfhist.nlambda;
+ buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
}
dd_bcast(dd, NITEM_DD_INIT_LOCAL_STATE*sizeof(int), buf);
t_atoms atoms;
atoms.nr = mtop->natoms;
- atoms.atom = NULL;
+ atoms.atom = nullptr;
make_reverse_ilist(mtop->intermolecular_ilist, &atoms,
- NULL, FALSE, FALSE, FALSE, TRUE, &ril_intermol);
+ nullptr, FALSE, FALSE, FALSE, TRUE, &ril_intermol);
}
snew(link, 1);
snew(link->index, ncg_mtop(mtop)+1);
link->nalloc_a = 0;
- link->a = NULL;
+ link->a = nullptr;
link->index[0] = 0;
cg_offset = 0;
* The constraints are discarded here.
*/
make_reverse_ilist(molt->ilist, &molt->atoms,
- NULL, FALSE, FALSE, FALSE, TRUE, &ril);
+ nullptr, FALSE, FALSE, FALSE, TRUE, &ril);
cgi_mb = &cginfo_mb[mb];
if (ePBC != epbcNONE)
{
- mk_mshift(NULL, graph, ePBC, box, x);
+ mk_mshift(nullptr, graph, ePBC, box, x);
shift_x(graph, box, x, xs);
/* By doing an extra mk_mshift the molecules that are broken
* will be made whole again. Such are the healing powers
* of GROMACS.
*/
- mk_mshift(NULL, graph, ePBC, box, xs);
+ mk_mshift(nullptr, graph, ePBC, box, xs);
}
else
{
if (vsite)
{
- construct_vsites(vsite, xs, 0.0, NULL,
+ construct_vsites(vsite, xs, 0.0, nullptr,
ffparams->iparams, molt->ilist,
- epbcNONE, TRUE, NULL, NULL);
+ epbcNONE, TRUE, nullptr, nullptr);
}
- calc_cgcm(NULL, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm);
+ calc_cgcm(nullptr, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm);
}
//! Returns whether \p molt has a virtual site
bExclRequired = inputrecExclForces(ir);
- vsite = init_vsite(mtop, NULL, TRUE);
+ vsite = init_vsite(mtop, nullptr, TRUE);
*r_2b = 0;
*r_mb = 0;
{
if (ir->ePBC != epbcNONE)
{
- mk_graph_ilist(NULL, molt->ilist, 0, molt->atoms.nr, FALSE, FALSE,
+ mk_graph_ilist(nullptr, molt->ilist, 0, molt->atoms.nr, FALSE, FALSE,
&graph);
}
for (mol = 0; mol < molb->nmol; mol++)
{
get_cgcm_mol(molt, &mtop->ffparams, ir->ePBC, &graph, box,
- have_vsite_molt(molt) ? vsite : NULL,
+ have_vsite_molt(molt) ? vsite : nullptr,
x+at_offset, xs, cg_cm);
bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
*/
#include "gmxpre.h"
+#include <cassert>
#include <cmath>
#include <cstdlib>
#include <cstring>
}
else
{
+ assert(sy);
ymin = sy[g][0][0];
ymax = sy[g][0][0];
for (s = 0; s < nsetspergraph; s++)
// better to be safe than sorry, so we check it with an assert.
// If we are in this comparison routine in the first place, neig2 should not be 0,
// so eigval2 should always be a valid pointer.
- GMX_RELEASE_ASSERT(eigval2 != NULL, "NULL pointer provided for eigval2");
+ GMX_RELEASE_ASSERT(eigval2 != nullptr, "NULL pointer provided for eigval2");
for (i = n; i < neig2; i++)
{
int noutvec, int *outvec, gmx_bool bSplit,
const gmx_output_env_t *oenv)
{
- FILE *xvgrout = NULL;
+ FILE *xvgrout = nullptr;
int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
- t_trxstatus *out = NULL;
+ t_trxstatus *out = nullptr;
t_trxstatus *status;
int noutvec_extr, imin, imax;
real *pmin, *pmax;
int *all_at;
matrix box;
rvec *xread, *x;
- real t, inp, **inprod = NULL;
+ real t, inp, **inprod = nullptr;
char str[STRLEN], str2[STRLEN], **ylabel, *c;
real fact;
- gmx_rmpbc_t gpbc = NULL;
+ gmx_rmpbc_t gpbc = nullptr;
snew(x, natoms);
/* calculate x: a fitted struture of the selected atoms */
if (bFit)
{
- reset_x(nfit, ifit, nat, NULL, xread, w_rls);
+ reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
do_fit(nat, w_rls, xref, xread);
}
for (i = 0; i < natoms; i++)
}
}
}
- write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
+ write_trx(out, natoms, index, atoms, 0, t, box, xread, nullptr, nullptr);
}
nframes++;
}
if (projfile)
{
- GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL if projfile is non-NULL");
+ GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL if projfile is non-NULL");
snew(ylabel, noutvec);
for (v = 0; v < noutvec; v++)
{
ylabel[v] = gmx_strdup(str);
}
sprintf(str, "projection on eigenvectors (%s)", proj_unit);
- write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
+ write_xvgr_graphs(projfile, noutvec, 1, str, nullptr, output_env_get_xvgr_tlabel(oenv),
(const char **)ylabel,
- nframes, inprod[noutvec], inprod, NULL,
+ nframes, inprod[noutvec], inprod, nullptr,
output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
}
{
t_atoms atoms;
rvec *x;
- real *b = NULL;
+ real *b = nullptr;
matrix box;
char *resnm, *atnm;
gmx_bool bPDB, b4D;
}
if ( ( b4D || bSplit ) && bPDB)
{
- GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL with 4D or split PDB output options");
+ GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL with 4D or split PDB output options");
out = gmx_ffopen(threedplotfile, "w");
fprintf(out, "HEADER %s\n", str);
}
else
{
- write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
+ write_sto_conf(threedplotfile, str, &atoms, x, nullptr, ePBC, box);
}
done_atom(&atoms);
}
snew(pmax, noutvec_extr);
if (extreme == 0)
{
- GMX_RELEASE_ASSERT(inprod != NULL, "inprod must be non-NULL");
+ GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
fprintf(stderr,
"%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
*eigvec[outvec[v]][i][d]/sqrtm[i]);
}
}
- write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
+ write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
}
close_trx(out);
}
write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
"black: total, red: x, green: y, blue: z",
"Atom number", (const char **)ylabel,
- natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
+ natoms, x, nullptr, y, 1, FALSE, FALSE, oenv);
fprintf(stderr, "\n");
}
y[g][i] = std::sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
}
}
- write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
+ write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", nullptr,
"Atom number", (const char **)ylabel,
- natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
+ natoms, x, y, nullptr, 1, TRUE, FALSE, oenv);
fprintf(stderr, "\n");
}
"of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
"have been set explicitly.[PAR]",
- "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
- "a single number for the overlap between the covariance matrices is",
- "generated. The formulas are::",
+ "When [TT]-v[tt] and [TT]-v2[tt] are given, a single number for the",
+ "overlap between the covariance matrices is generated. Note that the",
+ "eigenvalues are by default read from the timestamp field in the",
+ "eigenvector input files, but when [TT]-eig[tt], or [TT]-eig2[tt] are",
+ "given, the corresponding eigenvalues are used instead. The formulas are::",
"",
" difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
" normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
t_topology top;
int ePBC = -1;
- const t_atoms *atoms = NULL;
- rvec *xtop, *xref1, *xref2, *xrefp = NULL;
+ const t_atoms *atoms = nullptr;
+ rvec *xtop, *xref1, *xref2, *xrefp = nullptr;
gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
- int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
- rvec *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
+ int nvec1, nvec2, *eignr1 = nullptr, *eignr2 = nullptr;
+ rvec *xav1, *xav2, **eigvec1 = nullptr, **eigvec2 = nullptr;
matrix topbox;
real totmass, *sqrtm, *w_rls, t;
int natoms;
const char *indexfile;
int i, j, d;
int nout, *iout, noutvec, *outvec, nfit;
- int *index = NULL, *ifit = NULL;
+ int *index = nullptr, *ifit = nullptr;
const char *VecFile, *Vec2File, *topfile;
const char *EigFile, *Eig2File;
const char *CompFile, *RmsfFile, *ProjOnVecFile;
const char *OverlapFile, *InpMatFile;
gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
- real *eigval1 = NULL, *eigval2 = NULL;
+ real *eigval1 = nullptr, *eigval2 = nullptr;
int neig1, neig2;
double **xvgdata;
gmx_output_env_t *oenv;
t_filenm fnm[] = {
{ efTRN, "-v", "eigenvec", ffREAD },
{ efTRN, "-v2", "eigenvec2", ffOPTRD },
- { efTRX, "-f", NULL, ffOPTRD },
- { efTPS, NULL, NULL, ffOPTRD },
- { efNDX, NULL, NULL, ffOPTRD },
+ { efTRX, "-f", nullptr, ffOPTRD },
+ { efTPS, nullptr, nullptr, ffOPTRD },
+ { efNDX, nullptr, nullptr, ffOPTRD },
{ efXVG, "-eig", "eigenval", ffOPTRD },
{ efXVG, "-eig2", "eigenval2", ffOPTRD },
{ efXVG, "-comp", "eigcomp", ffOPTWR },
if (!parse_common_args(&argc, argv,
PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
- NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
+ NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
neig1 = DIM*natoms;
/* Overwrite eigenvalues from separate files if the user provides them */
- if (EigFile != NULL)
+ if (EigFile != nullptr)
{
int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
if (neig_tmp != neig1)
neig2 = 0;
}
- if (Eig2File != NULL)
+ if (Eig2File != nullptr)
{
neig2 = read_xvg(Eig2File, &xvgdata, &i);
srenew(eigval2, neig2);
{
bM = FALSE;
}
- if ((xref1 == NULL) && (bM || bTraj))
+ if ((xref1 == nullptr) && (bM || bTraj))
{
bTPS = TRUE;
}
- xtop = NULL;
+ xtop = nullptr;
nfit = 0;
- ifit = NULL;
- w_rls = NULL;
+ ifit = nullptr;
+ w_rls = nullptr;
if (!bTPS)
{
else
{
bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
- &top, &ePBC, &xtop, NULL, topbox, bM);
+ &top, &ePBC, &xtop, nullptr, topbox, bM);
atoms = &top.atoms;
gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
/* Fitting is only required for the projection */
if (bProj && bFit1)
{
- if (xref1 == NULL)
+ if (xref1 == nullptr)
{
printf("\nNote: the structure in %s should be the same\n"
" as the one used for the fit in g_covar\n", topfile);
}
snew(xrefp, atoms->nr);
- if (xref1 != NULL)
+ if (xref1 != nullptr)
{
/* Safety check between selected fit-group and reference structure read from the eigenvector file */
if (natoms != nfit)
{
copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
}
- reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
+ reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
}
}
gmx_rmpbc_done(gpbc);
{
printf("Select eigenvectors for output, end your selection with 0\n");
nout = -1;
- iout = NULL;
+ iout = nullptr;
do
{
if (bProj)
{
- project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
- bTop ? &top : NULL, ePBC, topbox,
+ project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr,
+ bTop ? &top : nullptr, ePBC, topbox,
ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
bFit1, xrefp, nfit, ifit, w_rls,
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
++ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
"the reference structure for the fit is written first with t=-1.",
"The average (or reference when [TT]-ref[tt] is used) structure is",
"written with t=0, the eigenvectors",
- "are written as frames with the eigenvector number as timestamp.",
+ "are written as frames with the eigenvector number and eigenvalue",
+ "as step number and timestamp, respectively.",
"[PAR]",
"The eigenvectors can be analyzed with [gmx-anaeig].",
"[PAR]",
{ "-pbc", FALSE, etBOOL, {&bPBC},
"Apply corrections for periodic boundary conditions" }
};
- FILE *out = NULL; /* initialization makes all compilers happy */
+ FILE *out = nullptr; /* initialization makes all compilers happy */
t_trxstatus *status;
t_topology top;
int ePBC;
matrix box, zerobox;
real *sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
real t, tstart, tend, **mat2;
- real xj, *w_rls = NULL;
+ real xj, *w_rls = nullptr;
real min, max, *axis;
int natoms, nat, nframes0, nframes, nlevels;
gmx_int64_t ndim, i, j, k, l;
t_rgb rlo, rmi, rhi;
real *eigenvectors;
gmx_output_env_t *oenv;
- gmx_rmpbc_t gpbc = NULL;
+ gmx_rmpbc_t gpbc = nullptr;
t_filenm fnm[] = {
- { efTRX, "-f", NULL, ffREAD },
- { efTPS, NULL, NULL, ffREAD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efXVG, NULL, "eigenval", ffWRITE },
+ { efTRX, "-f", nullptr, ffREAD },
+ { efTPS, nullptr, nullptr, ffREAD },
+ { efNDX, nullptr, nullptr, ffOPTRD },
+ { efXVG, nullptr, "eigenval", ffWRITE },
{ efTRN, "-v", "eigenvec", ffWRITE },
{ efSTO, "-av", "average.pdb", ffWRITE },
- { efLOG, NULL, "covar", ffWRITE },
+ { efLOG, nullptr, "covar", ffWRITE },
{ efDAT, "-ascii", "covar", ffOPTWR },
{ efXPM, "-xpm", "covar", ffOPTWR },
{ efXPM, "-xpma", "covara", ffOPTWR }
#define NFILE asize(fnm)
if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
xpmfile = opt2fn_null("-xpm", NFILE, fnm);
xpmafile = opt2fn_null("-xpma", NFILE, fnm);
- read_tps_conf(fitfile, &top, &ePBC, &xref, NULL, box, TRUE);
+ read_tps_conf(fitfile, &top, &ePBC, &xref, nullptr, box, TRUE);
atoms = &top.atoms;
if (bFit)
}
if (bFit)
{
- reset_x(nfit, ifit, atoms->nr, NULL, xref, w_rls);
+ reset_x(nfit, ifit, atoms->nr, nullptr, xref, w_rls);
}
snew(x, natoms);
}
if (bFit)
{
- reset_x(nfit, ifit, nat, NULL, xread, w_rls);
+ reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
do_fit(nat, w_rls, xref, xread);
}
for (i = 0; i < natoms; i++)
}
}
write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
- atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
+ atoms, xread, nullptr, epbcNONE, zerobox, natoms, index);
sfree(xread);
fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", static_cast<int>(ndim), static_cast<int>(ndim));
}
if (bFit)
{
- reset_x(nfit, ifit, nat, NULL, xread, w_rls);
+ reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
do_fit(nat, w_rls, xref, xread);
}
if (bRef)
int i, n;
gmx_bool found;
+ if (n_bfac > atoms->nres)
+ {
+ peratom = TRUE;
+ }
+
bfac_max = -1e10;
bfac_min = 1e10;
for (i = 0; (i < n_bfac); i++)
{
- if (bfac_nr[i] - 1 >= atoms->nres)
- {
- peratom = TRUE;
- }
/* if ((bfac_nr[i]-1<0) || (bfac_nr[i]-1>=atoms->nr))
gmx_fatal(FARGS,"Index of B-Factor %d is out of range: %d (%g)",
i+1,bfac_nr[i],bfac[i]); */
}
calc_triclinic_images(box, img + 1);
- write_sto_conf(fn, "Images", &atoms, img, NULL, ePBC, box);
+ write_sto_conf(fn, "Images", &atoms, img, nullptr, ePBC, box);
done_atom(&atoms);
sfree(img);
resind_prev = -1;
for (i = 0; i < isize; i++)
{
- resind = atoms->atom[index == NULL ? i : index[i]].resind;
+ resind = atoms->atom[index == nullptr ? i : index[i]].resind;
if (resind != resind_prev)
{
atoms->resinfo[resind].nr = resnr_start;
"from a file with with following format: first line states number of",
"entries in the file, next lines state an index",
"followed by a B-factor. The B-factors will be attached per residue",
- "unless an index is larger than the number of residues or unless the",
+ "unless the number of B-factors is larger than the number of the residues or unless the",
"[TT]-atom[tt] option is set. Obviously, any type of numeric data can",
"be added instead of B-factors. [TT]-legend[tt] will produce",
"a row of CA atoms with B-factors ranging from the minimum to the",
{ 0, 0, 0 }, targetvec =
{ 0, 0, 0 };
static const char *btype[] =
- { NULL, "triclinic", "cubic", "dodecahedron", "octahedron", NULL },
+ { nullptr, "triclinic", "cubic", "dodecahedron", "octahedron", nullptr },
*label = "A";
static rvec visbox =
{ 0, 0, 0 };
FILE *out;
const char *infile, *outfile;
int outftp, inftp, natom, i, j, n_bfac, itype, ntype;
- double *bfac = NULL, c6, c12;
- int *bfac_nr = NULL;
- t_topology *top = NULL;
+ double *bfac = nullptr, c6, c12;
+ int *bfac_nr = nullptr;
+ t_topology *top = nullptr;
char *grpname, *sgrpname, *agrpname;
int isize, ssize, asize;
int *index, *sindex, *aindex;
gmx_output_env_t *oenv;
t_filenm fnm[] =
{
- { efSTX, "-f", NULL, ffREAD },
- { efNDX, "-n", NULL, ffOPTRD },
- { efSTO, NULL, NULL, ffOPTWR },
+ { efSTX, "-f", nullptr, ffREAD },
+ { efNDX, "-n", nullptr, ffOPTRD },
+ { efSTO, nullptr, nullptr, ffOPTWR },
{ efPQR, "-mead", "mead", ffOPTWR },
{ efDAT, "-bf", "bfact", ffOPTRD }
};
bScale = bScale || bRho;
bCalcGeom = bCenter || bRotate || bOrient || bScale;
- GMX_RELEASE_ASSERT(btype[0] != NULL, "Option setting inconsistency; btype[0] is NULL");
+ GMX_RELEASE_ASSERT(btype[0] != nullptr, "Option setting inconsistency; btype[0] is NULL");
bCalcDiam = (btype[0][0] == 'c' || btype[0][0] == 'd' || btype[0][0] == 'o');
read_tps_conf(infile, top_tmp, &ePBC, &x, &v, box, FALSE);
t_atoms &atoms = top_tmp->atoms;
natom = atoms.nr;
- if (atoms.pdbinfo == NULL)
+ if (atoms.pdbinfo == nullptr)
{
snew(atoms.pdbinfo, atoms.nr);
}
+ atoms.havePdbInfo = TRUE;
+
if (fn2ftp(infile) == efPDB)
{
get_pdb_atomnumber(&atoms, aps);
if (bMead || bGrasp || bCONECT)
{
- top = read_top(infile, NULL);
+ top = read_top(infile, nullptr);
}
if (bMead || bGrasp)
else
{
ssize = atoms.nr;
- sindex = NULL;
+ sindex = nullptr;
}
diam = calc_geom(ssize, sindex, x, gc, rmin, rmax, bCalcDiam);
rvec_sub(rmax, rmin, size);
get_index(&atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &isize, &index, &grpnames);
/* Orient the principal axes along the coordinate axes */
- orient_princ(&atoms, isize, index, natom, x, bHaveV ? v : NULL, NULL);
+ orient_princ(&atoms, isize, index, natom, x, bHaveV ? v : nullptr, nullptr);
sfree(index);
sfree(grpnames);
}
else
{
ssize = atoms.nr;
- sindex = NULL;
+ sindex = nullptr;
}
printf("Translating %d atoms (out of %d) by %g %g %g nm\n", ssize, natom,
translation[XX], translation[YY], translation[ZZ]);
}
}
- if ((btype[0] != NULL) && (bSetSize || bDist || (btype[0][0] == 't' && bSetAng)))
+ if ((btype[0] != nullptr) && (bSetSize || bDist || (btype[0][0] == 't' && bSetAng)))
{
ePBC = epbcXYZ;
if (!(bSetSize || bDist))
}
else
{
- conect = NULL;
+ conect = nullptr;
}
if (bIndex)
}
else
{
- write_sto_conf_indexed(outfile, *top_tmp->name, &atoms, x, bHaveV ? v : NULL, ePBC, box, isize, index);
+ write_sto_conf_indexed(outfile, *top_tmp->name, &atoms, x, bHaveV ? v : nullptr, ePBC, box, isize, index);
}
}
else
{
if (resnr_start >= 0)
{
- renum_resnr(&atoms, atoms.nr, NULL, resnr_start);
+ renum_resnr(&atoms, atoms.nr, nullptr, resnr_start);
}
if ((outftp == efPDB) || (outftp == efPQR))
}
else
{
- write_sto_conf(outfile, *top_tmp->name, &atoms, x, bHaveV ? v : NULL, ePBC, box);
+ write_sto_conf(outfile, *top_tmp->name, &atoms, x, bHaveV ? v : nullptr, ePBC, box);
}
}
gmx_atomprop_destroy(aps);
- do_view(oenv, outfile, NULL);
+ do_view(oenv, outfile, nullptr);
return 0;
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
eigensolver(hess, ndim, begin-1, end-1, eigenvalues, eigenvectors);
/* And scale the output eigenvectors */
- if (bM && eigenvectors != NULL)
+ if (bM && eigenvectors != nullptr)
{
for (int i = 0; i < (end-begin+1); i++)
{
/* Cannot check symmetry since we only store half matrix */
/* divide elements hess[i][j] by sqrt(mas[i])*sqrt(mas[j]) when required */
- GMX_RELEASE_ASSERT(sparse_hessian != NULL, "NULL matrix pointer provided to nma_sparse_hessian");
+ GMX_RELEASE_ASSERT(sparse_hessian != nullptr, "NULL matrix pointer provided to nma_sparse_hessian");
if (bM)
{
sparse_eigensolver(sparse_hessian, neig, eigenvalues, eigenvectors, 10000000);
/* Scale output eigenvectors */
- if (bM && eigenvectors != NULL)
+ if (bM && eigenvectors != nullptr)
{
for (i = 0; i < neig; i++)
{
"which can be calculated with [gmx-mdrun].",
"The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
"The structure is written first with t=0. The eigenvectors",
- "are written as frames with the eigenvector number as timestamp.",
+ "are written as frames with the eigenvector number and eigenvalue",
+ "written as step number and timestamp, respectively.",
"The eigenvectors can be analyzed with [gmx-anaeig].",
"An ensemble of structures can be generated from the eigenvectors with",
"[gmx-nmens]. When mass weighting is used, the generated eigenvectors",
real value, omega, nu;
real factor_gmx_to_omega2;
real factor_omega_to_wavenumber;
- real *spectrum = NULL;
+ real *spectrum = nullptr;
real wfac;
gmx_output_env_t *oenv;
const char *qcleg[] = {
"Heat Capacity cV (J/mol K)",
"Enthalpy H (kJ/mol)"
};
- real * full_hessian = NULL;
- gmx_sparsematrix_t * sparse_hessian = NULL;
+ real * full_hessian = nullptr;
+ gmx_sparsematrix_t * sparse_hessian = nullptr;
t_filenm fnm[] = {
{ efMTX, "-f", "hessian", ffREAD },
- { efTPR, NULL, NULL, ffREAD },
+ { efTPR, nullptr, nullptr, ffREAD },
{ efXVG, "-of", "eigenfreq", ffWRITE },
{ efXVG, "-ol", "eigenval", ffWRITE },
{ efXVG, "-os", "spectrum", ffOPTWR },
#define NFILE asize(fnm)
if (!parse_common_args(&argc, argv, 0,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
{
return 0;
}
snew(top_x, tpx.natoms);
int natoms_tpx;
- read_tpx(ftp2fn(efTPR, NFILE, fnm), NULL, box, &natoms_tpx,
- top_x, NULL, &mtop);
+ read_tpx(ftp2fn(efTPR, NFILE, fnm), nullptr, box, &natoms_tpx,
+ top_x, nullptr, &mtop);
int nharm = 0;
if (bCons)
{
}
std::vector<size_t> atom_index = get_atom_index(&mtop);
- top = gmx_mtop_t_to_t_topology(&mtop);
+ top = gmx_mtop_t_to_t_topology(&mtop, true);
bM = TRUE;
int ndim = DIM*atom_index.size();
* If this is not valid we convert to full matrix storage,
* but warn the user that we might run out of memory...
*/
- if ((sparse_hessian != NULL) && (end == ndim))
+ if ((sparse_hessian != nullptr) && (end == ndim))
{
fprintf(stderr, "Cannot use sparse Hessian to calculate all eigenvectors.\n");
}
}
gmx_sparsematrix_destroy(sparse_hessian);
- sparse_hessian = NULL;
+ sparse_hessian = nullptr;
fprintf(stderr, "Converted sparse to full matrix storage.\n");
}
snew(eigenvalues, nrow);
- if (full_hessian != NULL)
+ if (full_hessian != nullptr)
{
/* Using full matrix storage */
eigenvectors = allocateEigenvectors(nrow, begin, end, false);
}
else
{
- qc = NULL;
+ qc = nullptr;
}
printf("Writing eigenfrequencies - negative eigenvalues will be set to zero.\n");
}
}
/* Spectrum ? */
- spec = NULL;
+ spec = nullptr;
if (opt2bSet("-os", NFILE, fnm) && (maxspec > 0))
{
snew(spectrum, maxspec);
nu = 1e-12*omega/(2*M_PI);
value = omega*factor_omega_to_wavenumber;
fprintf (out, "%6d %15g\n", i, value);
- if (NULL != spec)
+ if (nullptr != spec)
{
wfac = eigenvalues[i-begin]/(width*std::sqrt(2*M_PI));
for (j = 0; (j < maxspec); j++)
spectrum[j] += wfac*std::exp(-gmx::square(j-value)/(2*gmx::square(width)));
}
}
- if (NULL != qc)
+ if (nullptr != qc)
{
qcv = cv_corr(nu, T);
qu = u_corr(nu, T);
}
}
xvgrclose(out);
- if (NULL != spec)
+ if (nullptr != spec)
{
for (j = 0; (j < maxspec); j++)
{
}
xvgrclose(spec);
}
- if (NULL != qc)
+ if (nullptr != qc)
{
printf("Quantum corrections for harmonic degrees of freedom\n");
printf("Use appropriate -first and -last options to get reliable results.\n");
* will not be strictly orthogonal in plain cartesian scalar products.
*/
const real *eigenvectorPtr;
- if (full_hessian != NULL)
+ if (full_hessian != nullptr)
{
eigenvectorPtr = eigenvectors;
}
eigenvectorPtr = eigenvectors + (begin - 1)*atom_index.size();
}
write_eigenvectors(opt2fn("-v", NFILE, fnm), atom_index.size(), eigenvectorPtr, FALSE, begin, end,
- eWXR_NO, NULL, FALSE, top_x, bM, eigenvalues);
+ eWXR_NO, nullptr, FALSE, top_x, bM, eigenvalues);
return 0;
}
#include "gromacs/commandline/pargs.h"
#include "gromacs/fft/calcgrid.h"
#include "gromacs/fileio/checkpoint.h"
-#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/tpxio.h"
#include "gromacs/gmxana/gmx_ana.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/perf_est.h"
+#include "gromacs/mdrunutility/mdmodules.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
fp = fopen(fn_out, "r");
fprintf(stdout, "\n\n");
- while (fgets(buf, STRLEN-1, fp) != NULL)
+ while (fgets(buf, STRLEN-1, fp) != nullptr)
{
fprintf(stdout, "%s", buf);
}
iFound = eFoundDDStr; /* Skip some case statements */
}
- while (fgets(line, STRLEN, fp) != NULL)
+ while (fgets(line, STRLEN, fp) != nullptr)
{
/* Remove leading spaces */
ltrim(line);
/* Check for TERM and INT signals from user: */
- if (std::strstr(line, errSIG) != NULL)
+ if (std::strstr(line, errSIG) != nullptr)
{
fclose(fp);
cleandata(perfdata, test_nr);
/* Check whether cycle resetting worked */
if (presteps > 0 && !bFoundResetStr)
{
- if (std::strstr(line, matchstrcr) != NULL)
+ if (std::strstr(line, matchstrcr) != nullptr)
{
sprintf(dumstring, "step %s", "%" GMX_SCNd64);
sscanf(line, dumstring, &resetsteps);
if (gmx_fexist(errfile))
{
fp = fopen(errfile, "r");
- while (fgets(line, STRLEN, fp) != NULL)
+ while (fgets(line, STRLEN, fp) != nullptr)
{
if (str_starts(line, "Fatal error:") )
{
- if (fgets(line, STRLEN, fp) != NULL)
+ if (fgets(line, STRLEN, fp) != nullptr)
{
fprintf(stderr, "\nWARNING: An error occurred during this benchmark:\n"
"%s\n", line);
/* Get the commands we need to set up the runs from environment variables */
if (!bThreads)
{
- if ( (cp = getenv("MPIRUN")) != NULL)
+ if ( (cp = getenv("MPIRUN")) != nullptr)
{
*cmd_mpirun = gmx_strdup(cp);
}
*cmd_mpirun = gmx_strdup(empty_mpirun);
}
- if (*cmd_mdrun == NULL)
+ if (*cmd_mdrun == nullptr)
{
/* The use of MDRUN is deprecated, but made available in 5.1
for backward compatibility. It may be removed in a future
version. */
- if ( (cp = getenv("MDRUN" )) != NULL)
+ if ( (cp = getenv("MDRUN" )) != nullptr)
{
*cmd_mdrun = gmx_strdup(cp);
}
const char *cmd_mdrun,
gmx_bool bNeedGpuSupport)
{
- char *command = NULL;
+ char *command = nullptr;
char *cp;
char line[STRLEN];
FILE *fp;
while (!feof(fp) )
{
cp = fgets(line, STRLEN, fp);
- if (cp != NULL)
+ if (cp != nullptr)
{
if (str_starts(line, match_mdrun) )
{
const char *fn_best_tpr, /* tpr file with the best performance */
const char *fn_sim_tpr) /* name of tpr file to be launched */
{
- t_inputrec *ir;
- t_state state;
- gmx_mtop_t mtop;
- char buf[200];
+ t_inputrec *ir;
+ t_state state;
+ gmx_mtop_t mtop;
+ char buf[200];
- snew(ir, 1);
+ gmx::MDModules mdModules;
+ ir = mdModules.inputrec();
read_tpx_state(fn_best_tpr, ir, &state, &mtop);
/* Reset nsteps and init_step to the value of the input .tpr file */
fprintf(stdout, buf, ir->nsteps);
fflush(stdout);
write_tpx_state(fn_sim_tpr, ir, &state, &mtop);
-
- sfree(ir);
}
static gmx_bool can_scale_rvdw(int vdwtype)
}
fprintf(stdout, ".\n");
-
- snew(ir, 1);
+ gmx::MDModules mdModules;
+ ir = mdModules.inputrec();
read_tpx_state(fn_sim_tpr, ir, &state, &mtop);
/* Check if some kind of PME was chosen */
/* Scale the Fourier grid spacing */
ir->nkx = ir->nky = ir->nkz = 0;
- calc_grid(NULL, state.box, fourierspacing*fac, &ir->nkx, &ir->nky, &ir->nkz);
+ calc_grid(nullptr, state.box, fourierspacing*fac, &ir->nkx, &ir->nky, &ir->nkz);
/* Adjust other radii since various conditions need to be fulfilled */
if (eelPME == ir->coulombtype)
}
fflush(stdout);
fflush(fp);
- sfree(ir);
}
{
char numstring[STRLEN];
char newfilename[STRLEN];
- const char *fn = NULL;
+ const char *fn = nullptr;
int i;
const char *opt;
* constructing mdrun command lines */
{
int i, nr, k, ret, count = 0, totaltests;
- int *nPMEnodes = NULL;
- t_perf *pd = NULL;
+ int *nPMEnodes = nullptr;
+ t_perf *pd = nullptr;
int cmdline_length;
char *command, *cmd_stub;
char buf[STRLEN];
/* Loop over various numbers of PME nodes: */
for (i = 0; i < *pmeentries; i++)
{
- char *cmd_gpu_ids = NULL;
+ char *cmd_gpu_ids = nullptr;
pd = &perfdata[k][i];
gmx_bool bFree; /* Is a free energy simulation requested? */
gmx_bool bNM; /* Is a normal mode analysis requested? */
gmx_bool bSwap; /* Is water/ion position swapping requested? */
- t_inputrec ir;
+ t_inputrec *ir;
t_state state;
gmx_mtop_t mtop;
/* Check tpr file for options that trigger extra output files */
- read_tpx_state(opt2fn("-s", nfile, fnm), &ir, &state, &mtop);
- bFree = (efepNO != ir.efep );
- bNM = (eiNM == ir.eI );
- bSwap = (eswapNO != ir.eSwapCoords);
- bTpi = EI_TPI(ir.eI);
+ gmx::MDModules mdModules;
+ ir = mdModules.inputrec();
+ read_tpx_state(opt2fn("-s", nfile, fnm), ir, &state, &mtop);
+ bFree = (efepNO != ir->efep );
+ bNM = (eiNM == ir->eI );
+ bSwap = (eswapNO != ir->eSwapCoords);
+ bTpi = EI_TPI(ir->eI);
/* Set these output files on the tuning command-line */
- if (ir.bPull)
+ if (ir->bPull)
{
setopt("-pf", nfile, fnm);
setopt("-px", nfile, fnm);
setopt("-swap", nfile, fnm);
}
- *rcoulomb = ir.rcoulomb;
+ *rcoulomb = ir->rcoulomb;
/* Return the estimate for the number of PME nodes */
- return pme_load_estimate(&mtop, &ir, state.box);
+ float npme = pme_load_estimate(&mtop, ir, state.box);
+ return npme;
}
real rcoulomb = -1.0; /* Coulomb radius as set in .tpr file */
gmx_bool bScaleRvdw = TRUE;
gmx_int64_t bench_nsteps = BENCHSTEPS;
- gmx_int64_t new_sim_nsteps = -1; /* -1 indicates: not set by the user */
- gmx_int64_t cpt_steps = 0; /* Step counter in .cpt input file */
- int presteps = 100; /* Do a full cycle reset after presteps steps */
+ gmx_int64_t new_sim_nsteps = -1; /* -1 indicates: not set by the user */
+ gmx_int64_t cpt_steps = 0; /* Step counter in .cpt input file */
+ int presteps = 1500; /* Do a full cycle reset after presteps steps */
gmx_bool bOverwrite = FALSE, bKeepTPR;
gmx_bool bLaunch = FALSE;
- char *ExtraArgs = NULL;
- char **tpr_names = NULL;
- const char *simulation_tpr = NULL;
- char *deffnm = NULL;
+ char *ExtraArgs = nullptr;
+ char **tpr_names = nullptr;
+ const char *simulation_tpr = nullptr;
+ char *deffnm = nullptr;
int best_npme, best_tpr;
int sim_part = 1; /* For benchmarks with checkpoint files */
char bbuf[STRLEN];
/* Default program names if nothing else is found */
- char *cmd_mpirun = NULL, *cmd_mdrun = NULL;
+ char *cmd_mpirun = nullptr, *cmd_mdrun = nullptr;
char *cmd_args_bench, *cmd_args_launch;
- char *cmd_np = NULL;
+ char *cmd_np = nullptr;
/* IDs of GPUs that are eligible for computation */
- char *eligible_gpu_ids = NULL;
- t_eligible_gpu_ids *gpu_ids = NULL;
+ char *eligible_gpu_ids = nullptr;
+ t_eligible_gpu_ids *gpu_ids = nullptr;
- t_perf **perfdata = NULL;
+ t_perf **perfdata = nullptr;
t_inputinfo *info;
int i;
FILE *fp;
{ efLOG, "-err", "bencherr", ffWRITE },
{ efTPR, "-so", "tuned", ffWRITE },
/* mdrun: */
- { efTPR, NULL, NULL, ffREAD },
- { efTRN, "-o", NULL, ffWRITE },
- { efCOMPRESSED, "-x", NULL, ffOPTWR },
- { efCPT, "-cpi", NULL, ffOPTRD },
- { efCPT, "-cpo", NULL, ffOPTWR },
+ { efTPR, nullptr, nullptr, ffREAD },
+ { efTRN, "-o", nullptr, ffWRITE },
+ { efCOMPRESSED, "-x", nullptr, ffOPTWR },
+ { efCPT, "-cpi", nullptr, ffOPTRD },
+ { efCPT, "-cpo", nullptr, ffOPTWR },
{ efSTO, "-c", "confout", ffWRITE },
{ efEDR, "-e", "ener", ffWRITE },
{ efLOG, "-g", "md", ffWRITE },
int nthreads = 1;
const char *procstring[] =
- { NULL, "np", "n", "none", NULL };
+ { nullptr, "np", "n", "none", nullptr };
const char *npmevalues_opt[] =
- { NULL, "auto", "all", "subset", NULL };
+ { nullptr, "auto", "all", "subset", nullptr };
gmx_bool bAppendFiles = TRUE;
gmx_bool bKeepAndNumCPT = FALSE;
gmx_bool bBenchmark = TRUE;
gmx_bool bCheck = TRUE;
- gmx_output_env_t *oenv = NULL;
+ gmx_output_env_t *oenv = nullptr;
t_pargs pa[] = {
/***********************/
if (!parse_common_args(&argc, argv, PCA_NOEXIT_ON_ARGS,
NFILE, fnm, asize(pa), pa, asize(desc), desc,
- 0, NULL, &oenv))
+ 0, nullptr, &oenv))
{
return 0;
}
// procstring[0]Â is used inside two different conditionals further down
- GMX_RELEASE_ASSERT(procstring[0] != NULL, "Options inconsistency; procstring[0]Â is NULL");
+ GMX_RELEASE_ASSERT(procstring[0] != nullptr, "Options inconsistency; procstring[0]Â is NULL");
/* Store the remaining unparsed command line entries in a string which
* is then attached to the mdrun command line */
get_program_paths(bThreads, &cmd_mpirun, &cmd_mdrun);
if (bBenchmark && repeats > 0)
{
- check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun, NULL != eligible_gpu_ids);
+ check_mdrun_works(bThreads, cmd_mpirun, cmd_np, cmd_mdrun, nullptr != eligible_gpu_ids);
}
/* Print some header info to file */
snew(perfdata, ntprs);
if (bBenchmark)
{
- GMX_RELEASE_ASSERT(npmevalues_opt[0] != NULL, "Options inconsistency; npmevalues_opt[0] is NULL");
+ GMX_RELEASE_ASSERT(npmevalues_opt[0] != nullptr, "Options inconsistency; npmevalues_opt[0] is NULL");
do_the_tests(fp, tpr_names, maxPMEnodes, minPMEnodes, npme_fixed, npmevalues_opt[0], perfdata, &pmeentries,
repeats, nnodes, ntprs, bThreads, cmd_mpirun, cmd_np, cmd_mdrun,
cmd_args_bench, fnm, NFILE, presteps, cpt_steps, bCheck, gpu_ids);
#include "gromacs/ewald/pme.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/mtxio.h"
-#include "gromacs/gmxlib/md_logging.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/imd/imd.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdebin.h"
#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/mdsetup.h"
#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/shellfc.h"
#include "gromacs/mdlib/sim_util.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/timing/walltime_accounting.h"
#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/smalloc.h"
//! Utility structure for manipulating states during EM
typedef struct {
//! Copy of the global state
- t_state s;
+ t_state s;
//! Force array
- rvec *f;
+ PaddedRVecVector f;
//! Potential energy
- real epot;
+ real epot;
//! Norm of the force
- real fnorm;
+ real fnorm;
//! Maximum force
- real fmax;
+ real fmax;
//! Direction
- int a_fmax;
+ int a_fmax;
} em_state_t;
-//! Initiate em_state_t structure and return pointer to it
-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;
-}
-
//! Print the EM starting conditions
static void print_em_start(FILE *fplog,
t_commrec *cr,
//! Print message about convergence of the EM
static void print_converged(FILE *fp, const char *alg, real ftol,
gmx_int64_t count, gmx_bool bDone, gmx_int64_t nsteps,
- real epot, real fmax, int nfmax, real fnorm)
+ const em_state_t *ems, double sqrtNumAtoms)
{
char buf[STEPSTRSIZE];
}
#if 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", ems->epot);
+ fprintf(fp, "Maximum force = %21.14e on atom %d\n", ems->fmax, ems->a_fmax + 1);
+ fprintf(fp, "Norm of force = %21.14e\n", ems->fnorm/sqrtNumAtoms);
#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", ems->epot);
+ fprintf(fp, "Maximum force = %14.7e on atom %d\n", ems->fmax, ems->a_fmax + 1);
+ fprintf(fp, "Norm of force = %14.7e\n", ems->fnorm/sqrtNumAtoms);
#endif
}
//! Compute the norm and max of the force array in parallel
static void get_f_norm_max(t_commrec *cr,
- t_grpopts *opts, t_mdatoms *mdatoms, rvec *f,
+ t_grpopts *opts, t_mdatoms *mdatoms, const rvec *f,
real *fnorm, real *fmax, int *a_fmax)
{
double fnorm2, *sum;
t_grpopts *opts, t_mdatoms *mdatoms,
em_state_t *ems)
{
- get_f_norm_max(cr, opts, mdatoms, ems->f, &ems->fnorm, &ems->fmax, &ems->a_fmax);
+ get_f_norm_max(cr, opts, mdatoms, as_rvec_array(ems->f.data()),
+ &ems->fnorm, &ems->fmax, &ems->a_fmax);
}
//! Initialize the energy minimization
t_commrec *cr, t_inputrec *ir,
t_state *state_global, gmx_mtop_t *top_global,
em_state_t *ems, gmx_localtop_t **top,
- rvec **f,
t_nrnb *nrnb, rvec mu_tot,
t_forcerec *fr, gmx_enerdata_t **enerd,
t_graph **graph, t_mdatoms *mdatoms, gmx_global_stat_t *gstat,
- gmx_vsite_t *vsite, gmx_constr_t constr,
+ gmx_vsite_t *vsite, gmx_constr_t constr, gmx_shellfc_t **shellfc,
int nfile, const t_filenm fnm[],
gmx_mdoutf_t *outf, t_mdebin **mdebin,
int imdport, unsigned long gmx_unused Flags,
gmx_wallcycle_t wcycle)
{
- int i;
real dvdl_constr;
if (fplog)
state_global->ngtc = 0;
/* Initialize lambda variables */
- initialize_lambdas(fplog, ir, &(state_global->fep_state), state_global->lambda, NULL);
+ initialize_lambdas(fplog, ir, &(state_global->fep_state), &state_global->lambda, nullptr);
init_nrnb(nrnb);
/* Interactive molecular dynamics */
- init_IMD(ir, cr, top_global, fplog, 1, state_global->x,
- nfile, fnm, NULL, imdport, Flags);
+ init_IMD(ir, cr, top_global, fplog, 1, as_rvec_array(state_global->x.data()),
+ nfile, fnm, nullptr, imdport, Flags);
+
+ if (ir->eI == eiNM)
+ {
+ GMX_ASSERT(shellfc != NULL, "With NM we always support shells");
+
+ *shellfc = init_shell_flexcon(stdout,
+ top_global,
+ n_flexible_constraints(constr),
+ ir->nstcalcenergy,
+ DOMAINDECOMP(cr));
+ }
+ else
+ {
+ GMX_ASSERT(EI_ENERGY_MINIMIZATION(ir->eI), "This else currently only handles energy minimizers, consider if your algorithm needs shell/flexible-constraint support");
+
+ /* With energy minimization, shells and flexible constraints are
+ * automatically minimized when treated like normal DOFS.
+ */
+ if (shellfc != nullptr)
+ {
+ *shellfc = nullptr;
+ }
+ }
if (DOMAINDECOMP(cr))
{
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,
&ems->s, &ems->f, mdatoms, *top,
fr, vsite, constr,
- nrnb, NULL, FALSE);
+ nrnb, nullptr, FALSE);
dd_store_state(cr->dd, &ems->s);
- *graph = NULL;
+ *graph = nullptr;
}
else
{
- snew(*f, top_global->natoms);
-
/* Just copy the state */
ems->s = *state_global;
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- snew(ems->s.x, ems->s.nalloc + 1);
- snew(ems->f, ems->s.nalloc+1);
- snew(ems->s.v, ems->s.nalloc+1);
- for (i = 0; i < state_global->natoms; i++)
- {
- copy_rvec(state_global->x[i], ems->s.x[i]);
- }
- copy_mat(state_global->box, ems->s.box);
+ ems->s.x.resize(ems->s.natoms + 1);
+ ems->f.resize(ems->s.natoms + 1);
- *top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
-
- setup_bonded_threading(fr, &(*top)->idef);
-
- if (ir->ePBC != epbcNONE && !fr->bMolPBC)
- {
- *graph = mk_graph(fplog, &((*top)->idef), 0, top_global->natoms, FALSE, FALSE);
- }
- else
- {
- *graph = NULL;
- }
-
- atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
+ snew(*top, 1);
+ mdAlgorithmsSetupAtomData(cr, ir, top_global, *top, fr,
+ graph, mdatoms,
+ vsite, shellfc ? *shellfc : nullptr);
- update_mdatoms(mdatoms, state_global->lambda[efptFEP]);
-
if (vsite)
{
set_vsite_top(vsite, *top, mdatoms, cr);
}
}
+ update_mdatoms(mdatoms, state_global->lambda[efptMASS]);
+
if (constr)
{
if (ir->eConstrAlg == econtSHAKE &&
{
/* Constrain the starting coordinates */
dvdl_constr = 0;
- constrain(PAR(cr) ? NULL : fplog, TRUE, TRUE, constr, &(*top)->idef,
+ constrain(PAR(cr) ? nullptr : fplog, TRUE, TRUE, constr, &(*top)->idef,
ir, cr, -1, 0, 1.0, mdatoms,
- ems->s.x, ems->s.x, NULL, fr->bMolPBC, ems->s.box,
+ as_rvec_array(ems->s.x.data()),
+ as_rvec_array(ems->s.x.data()),
+ nullptr,
+ fr->bMolPBC, ems->s.box,
ems->s.lambda[efptFEP], &dvdl_constr,
- NULL, NULL, nrnb, econqCoord);
+ nullptr, nullptr, nrnb, econqCoord);
}
}
}
else
{
- *gstat = NULL;
+ *gstat = nullptr;
}
- *outf = init_mdoutf(fplog, nfile, fnm, 0, cr, ir, top_global, NULL, wcycle);
+ *outf = init_mdoutf(fplog, nfile, fnm, 0, cr, ir, top_global, nullptr, wcycle);
snew(*enerd, 1);
init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
*enerd);
- if (mdebin != NULL)
+ if (mdebin != nullptr)
{
/* Init bin for energy stuff */
- *mdebin = init_mdebin(mdoutf_get_fp_ene(*outf), top_global, ir, NULL);
+ *mdebin = init_mdebin(mdoutf_get_fp_ene(*outf), top_global, ir, nullptr);
}
clear_rvec(mu_tot);
}
//! Finalize the minimization
-static void finish_em(t_commrec *cr, gmx_mdoutf_t outf,
+static void finish_em(t_commrec *cr, gmx_mdoutf_t outf, const t_inputrec *ir,
gmx_walltime_accounting_t walltime_accounting,
gmx_wallcycle_t wcycle)
{
gmx_pme_send_finish(cr);
}
- done_mdoutf(outf);
+ done_mdoutf(outf, ir);
em_time_end(walltime_accounting, wcycle);
}
//! Swap two different EM states during minimization
-static void swap_em_state(em_state_t *ems1, em_state_t *ems2)
+static void swap_em_state(em_state_t **ems1, em_state_t **ems2)
{
- em_state_t tmp;
+ em_state_t *tmp;
tmp = *ems1;
*ems1 = *ems2;
*ems2 = tmp;
}
-//! Copy coordinate from an EM state to a "normal" state structure
-static void copy_em_coords(em_state_t *ems, t_state *state)
-{
- int i;
-
- for (i = 0; (i < state->natoms); i++)
- {
- copy_rvec(ems->s.x[i], state->x[i]);
- }
-}
-
//! Save the EM trajectory
static void write_em_traj(FILE *fplog, t_commrec *cr,
gmx_mdoutf_t outf,
gmx_mtop_t *top_global,
t_inputrec *ir, gmx_int64_t step,
em_state_t *state,
- t_state *state_global)
+ t_state *state_global,
+ energyhistory_t *energyHistory)
{
- int mdof_flags;
- gmx_bool bIMDout = FALSE;
-
+ int mdof_flags = 0;
- /* Shall we do IMD output? */
- if (ir->bIMD)
- {
- bIMDout = do_per_step(step, IMD_get_step(ir->imd->setup));
- }
-
- if ((bX || bF || bIMDout || confout != NULL) && !DOMAINDECOMP(cr))
- {
- copy_em_coords(state, state_global);
- }
-
- mdof_flags = 0;
if (bX)
{
mdof_flags |= MDOF_X;
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
top_global, step, (double)step,
- &state->s, state_global, state->f);
+ &state->s, state_global, energyHistory,
+ &state->f);
- if (confout != NULL && MASTER(cr))
+ if (confout != nullptr && MASTER(cr))
{
if (ir->ePBC != epbcNONE && !ir->bPeriodicMols && DOMAINDECOMP(cr))
{
/* Make molecules whole only for confout writing */
do_pbc_mtop(fplog, ir->ePBC, state_global->box, top_global,
- state_global->x);
+ as_rvec_array(state_global->x.data()));
}
write_sto_conf_mtop(confout,
*top_global->name, top_global,
- state_global->x, NULL, ir->ePBC, state_global->box);
+ as_rvec_array(state_global->x.data()), nullptr, ir->ePBC, state_global->box);
}
}
// \returns true when the step succeeded, false when a constraint error occurred
static bool do_em_step(t_commrec *cr, t_inputrec *ir, t_mdatoms *md,
gmx_bool bMolPBC,
- em_state_t *ems1, real a, rvec *f, em_state_t *ems2,
+ em_state_t *ems1, real a, const PaddedRVecVector *force,
+ em_state_t *ems2,
gmx_constr_t constr, gmx_localtop_t *top,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_int64_t count)
s2->flags = s1->flags;
- if (s2->nalloc != s1->nalloc)
+ if (s2->natoms != s1->natoms)
{
- s2->nalloc = s1->nalloc;
+ s2->natoms = s1->natoms;
/* We need to allocate one element extra, since we might use
* (unaligned) 4-wide SIMD loads to access rvec entries.
*/
- srenew(s2->x, s1->nalloc + 1);
- srenew(ems2->f, s1->nalloc);
+ s2->x.resize(s1->natoms + 1);
+ ems2->f.resize(s1->natoms + 1);
if (s2->flags & (1<<estCGP))
{
- srenew(s2->cg_p, s1->nalloc + 1);
+ s2->cg_p.resize(s1->natoms + 1);
}
}
+ if (DOMAINDECOMP(cr) && s2->cg_gl.size() != s1->cg_gl.size())
+ {
+ s2->cg_gl.resize(s1->cg_gl.size());
+ }
s2->natoms = s1->natoms;
copy_mat(s1->box, s2->box);
/* Copy free energy state */
- for (int i = 0; i < efptNR; i++)
- {
- s2->lambda[i] = s1->lambda[i];
- }
+ s2->lambda = s1->lambda;
copy_mat(s1->box, s2->box);
start = 0;
nthreads = gmx_omp_nthreads_get(emntUpdate);
#pragma omp parallel num_threads(nthreads)
{
- rvec *x1 = s1->x;
- rvec *x2 = s2->x;
+ const rvec *x1 = as_rvec_array(s1->x.data());
+ rvec *x2 = as_rvec_array(s2->x.data());
+ const rvec *f = as_rvec_array(force->data());
- int gf = 0;
+ int gf = 0;
#pragma omp for schedule(static) nowait
for (int i = start; i < end; i++)
{
if (s2->flags & (1<<estCGP))
{
/* Copy the CG p vector */
- rvec *p1 = s1->cg_p;
- rvec *p2 = s2->cg_p;
+ const rvec *p1 = as_rvec_array(s1->cg_p.data());
+ rvec *p2 = as_rvec_array(s2->cg_p.data());
#pragma omp for schedule(static) nowait
for (int i = start; i < end; i++)
{
if (DOMAINDECOMP(cr))
{
s2->ddp_count = s1->ddp_count;
- if (s2->cg_gl_nalloc < s1->cg_gl_nalloc)
- {
-#pragma omp barrier
- s2->cg_gl_nalloc = s1->cg_gl_nalloc;
- try
- {
- /* We need to allocate one element extra, since we might use
- * (unaligned) 4-wide SIMD loads to access rvec entries.
- */
- srenew(s2->cg_gl, s2->cg_gl_nalloc + 1);
- }
- GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
-#pragma omp barrier
- }
- s2->ncg_gl = s1->ncg_gl;
+
+ /* OpenMP does not supported unsigned loop variables */
#pragma omp for schedule(static) nowait
- for (int i = 0; i < s2->ncg_gl; i++)
+ for (int i = 0; i < static_cast<int>(s2->cg_gl.size()); i++)
{
s2->cg_gl[i] = s1->cg_gl[i];
}
wallcycle_start(wcycle, ewcCONSTR);
dvdl_constr = 0;
validStep =
- constrain(NULL, TRUE, TRUE, constr, &top->idef,
+ constrain(nullptr, TRUE, TRUE, constr, &top->idef,
ir, cr, count, 0, 1.0, md,
- s1->x, s2->x, NULL, bMolPBC, s2->box,
+ as_rvec_array(s1->x.data()), as_rvec_array(s2->x.data()),
+ nullptr, bMolPBC, s2->box,
s2->lambda[efptBONDED], &dvdl_constr,
- NULL, NULL, nrnb, econqCoord);
+ nullptr, nullptr, nrnb, econqCoord);
wallcycle_stop(wcycle, ewcCONSTR);
// We should move this check to the different minimizers
{
/* Repartition the domain decomposition */
dd_partition_system(fplog, step, cr, FALSE, 1,
- NULL, top_global, ir,
+ nullptr, top_global, ir,
&ems->s, &ems->f,
mdatoms, top, fr, vsite, constr,
nrnb, wcycle, FALSE);
if (vsite)
{
- construct_vsites(vsite, ems->s.x, 1, NULL,
+ construct_vsites(vsite, as_rvec_array(ems->s.x.data()), 1, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, ems->s.box);
}
*/
do_force(fplog, cr, inputrec,
count, nrnb, wcycle, top, &top_global->groups,
- ems->s.box, ems->s.x, &ems->s.hist,
- ems->f, force_vir, mdatoms, enerd, fcd,
- ems->s.lambda, graph, fr, vsite, mu_tot, t, NULL, NULL, TRUE,
+ ems->s.box, &ems->s.x, &ems->s.hist,
+ &ems->f, force_vir, mdatoms, enerd, fcd,
+ &ems->s.lambda, graph, fr, vsite, mu_tot, t, nullptr, TRUE,
GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES |
GMX_FORCE_VIRIAL | GMX_FORCE_ENERGY |
(bNS ? GMX_FORCE_NS : 0));
wallcycle_start(wcycle, ewcMoveE);
global_stat(gstat, cr, enerd, force_vir, shake_vir, mu_tot,
- inputrec, NULL, NULL, NULL, 1, &terminate,
- NULL, FALSE,
+ inputrec, nullptr, nullptr, nullptr, 1, &terminate,
+ nullptr, FALSE,
CGLO_ENERGY |
CGLO_PRESSURE |
CGLO_CONSTRAINT);
/* Project out the constraint components of the force */
wallcycle_start(wcycle, ewcCONSTR);
dvdl_constr = 0;
- constrain(NULL, FALSE, FALSE, constr, &top->idef,
+ rvec *f_rvec = as_rvec_array(ems->f.data());
+ constrain(nullptr, FALSE, FALSE, constr, &top->idef,
inputrec, cr, count, 0, 1.0, mdatoms,
- ems->s.x, ems->f, ems->f, fr->bMolPBC, ems->s.box,
+ as_rvec_array(ems->s.x.data()), f_rvec, f_rvec,
+ fr->bMolPBC, ems->s.box,
ems->s.lambda[efptBONDED], &dvdl_constr,
- NULL, &shake_vir, nrnb, econqForceDispl);
+ nullptr, &shake_vir, nrnb, econqForceDispl);
enerd->term[F_DVDL_CONSTR] += dvdl_constr;
m_add(force_vir, shake_vir, vir);
wallcycle_stop(wcycle, ewcCONSTR);
enerd->term[F_PRES] =
calc_pres(fr->ePBC, inputrec->nwall, ems->s.box, ekin, vir, pres);
- sum_dhdl(enerd, ems->s.lambda, inputrec->fepvals);
+ sum_dhdl(enerd, &ems->s.lambda, inputrec->fepvals);
if (EI_ENERGY_MINIMIZATION(inputrec->eI))
{
gmx_mtop_t *top_global,
em_state_t *s_min, em_state_t *s_b)
{
- rvec *fm, *fb, *fmg;
t_block *cgs_gl;
int ncg, *cg_gl, *index, c, cg, i, a0, a1, a, gf, m;
double partsum;
fprintf(debug, "Doing reorder_partsum\n");
}
- fm = s_min->f;
- fb = s_b->f;
+ const rvec *fm = as_rvec_array(s_min->f.data());
+ const rvec *fb = as_rvec_array(s_b->f.data());
cgs_gl = dd_charge_groups_global(cr->dd);
index = cgs_gl->index;
* This conflicts with the spirit of domain decomposition,
* but to fully optimize this a much more complicated algorithm is required.
*/
+ rvec *fmg;
snew(fmg, top_global->natoms);
- ncg = s_min->s.ncg_gl;
- cg_gl = s_min->s.cg_gl;
+ ncg = s_min->s.cg_gl.size();
+ cg_gl = s_min->s.cg_gl.data();
i = 0;
for (c = 0; c < ncg; c++)
{
gmx_sum(top_global->natoms*3, fmg[0], cr);
/* Now we will determine the part of the sum for the cgs in state s_b */
- ncg = s_b->s.ncg_gl;
- cg_gl = s_b->s.cg_gl;
+ ncg = s_b->s.cg_gl.size();
+ cg_gl = s_b->s.cg_gl.data();
partsum = 0;
i = 0;
gf = 0;
gmx_mtop_t *top_global,
em_state_t *s_min, em_state_t *s_b)
{
- rvec *fm, *fb;
double sum;
- int gf, i, m;
/* This is just the classical Polak-Ribiere calculation of beta;
* it looks a bit complicated since we take freeze groups into account,
(s_min->s.ddp_count == cr->dd->ddp_count &&
s_b->s.ddp_count == cr->dd->ddp_count))
{
- fm = s_min->f;
- fb = s_b->f;
- sum = 0;
- gf = 0;
+ const rvec *fm = as_rvec_array(s_min->f.data());
+ const rvec *fb = as_rvec_array(s_b->f.data());
+ sum = 0;
+ int gf = 0;
/* This part of code can be incorrect with DD,
* since the atom ordering in s_b and s_min might differ.
*/
- for (i = 0; i < mdatoms->homenr; i++)
+ for (int i = 0; i < mdatoms->homenr; i++)
{
if (mdatoms->cFREEZE)
{
gf = mdatoms->cFREEZE[i];
}
- for (m = 0; m < DIM; m++)
+ for (int m = 0; m < DIM; m++)
{
if (!opts->nFreeze[gf][m])
{
{
/*! \brief Do conjugate gradients minimization
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+ \copydoc integrator_t(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t *oenv, gmx_bool bVerbose,
int nstglobalcomm,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
*/
-double do_cg(FILE *fplog, t_commrec *cr,
+double do_cg(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t gmx_unused *oenv, gmx_bool bVerbose,
int gmx_unused nstglobalcomm,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
+ energyhistory_t *energyHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_edsam_t gmx_unused ed,
{
const char *CG = "Polak-Ribiere Conjugate Gradients";
- em_state_t *s_min, *s_a, *s_b, *s_c;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
- rvec *p, *sf;
- double gpa, gpb, gpc, tmp, minstep;
- real fnormn;
+ double tmp, minstep;
real stepsize;
real a, b, c, beta = 0.0;
real epot_repl = 0;
tensor vir, pres;
int number_steps, neval = 0, nstcg = inputrec->nstcgsteep;
gmx_mdoutf_t outf;
- int i, m, gf, step, nminstep;
+ int m, step, nminstep;
step = 0;
- s_min = init_em_state();
- s_a = init_em_state();
- s_b = init_em_state();
- s_c = init_em_state();
+ /* Create 4 states on the stack and extract pointers that we will swap */
+ em_state_t s0 {}, s1 {}, s2 {}, s3 {};
+ em_state_t *s_min = &s0;
+ em_state_t *s_a = &s1;
+ em_state_t *s_b = &s2;
+ em_state_t *s_c = &s3;
/* Init em and store the local state in s_min */
init_em(fplog, CG, cr, inputrec,
- state_global, top_global, s_min, &top, &f,
- nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
+ state_global, top_global, s_min, &top,
+ nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
+ vsite, constr, nullptr,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
/* Print to log file */
/* Copy stuff to the energy bin for easy printing etc. */
upd_mdebin(mdebin, FALSE, FALSE, (double)step,
mdatoms->tmass, enerd, &s_min->s, inputrec->fepvals, inputrec->expandedvals, s_min->s.box,
- NULL, NULL, vir, pres, NULL, mu_tot, constr);
+ nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
print_ebin_header(fplog, step, step);
print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
*/
/* Calculate the new direction in p, and the gradient in this direction, gpa */
- p = s_min->s.cg_p;
- sf = s_min->f;
- gpa = 0;
- gf = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ rvec *pm = as_rvec_array(s_min->s.cg_p.data());
+ const rvec *sfm = as_rvec_array(s_min->f.data());
+ double gpa = 0;
+ int gf = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
if (mdatoms->cFREEZE)
{
{
if (!inputrec->opts.nFreeze[gf][m])
{
- p[i][m] = sf[i][m] + beta*p[i][m];
- gpa -= p[i][m]*sf[i][m];
+ pm[i][m] = sfm[i][m] + beta*pm[i][m];
+ gpa -= pm[i][m]*sfm[i][m];
/* f is negative gradient, thus the sign */
}
else
{
- p[i][m] = 0;
+ pm[i][m] = 0;
}
}
}
}
/* Calculate the norm of the search vector */
- get_f_norm_max(cr, &(inputrec->opts), mdatoms, p, &pnorm, NULL, NULL);
+ get_f_norm_max(cr, &(inputrec->opts), mdatoms, pm, &pnorm, nullptr, nullptr);
/* Just in case stepsize reaches zero due to numerical precision... */
if (stepsize <= 0)
* relative change in coordinate is smaller than precision
*/
minstep = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
{
tmp = 1.0;
}
- tmp = p[i][m]/tmp;
+ tmp = pm[i][m]/tmp;
minstep += tmp*tmp;
}
}
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,
+ write_em_traj(fplog, cr, outf, do_x, do_f, nullptr,
top_global, inputrec, step,
- s_min, state_global);
+ s_min, state_global, energyHistory);
/* Take a step downhill.
* In theory, we should minimize the function along this direction.
}
/* Take a trial step (new coords in s_c) */
- do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, c, s_min->s.cg_p, s_c,
+ do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, c, &s_min->s.cg_p, s_c,
constr, top, nrnb, wcycle, -1);
neval++;
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 = 0; i < mdatoms->homenr; i++)
+ const rvec *pc = as_rvec_array(s_c->s.cg_p.data());
+ const rvec *sfc = as_rvec_array(s_c->f.data());
+ double gpc = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
- gpc -= p[i][m]*sf[i][m]; /* f is negative gradient, thus the sign */
+ gpc -= pc[i][m]*sfc[i][m]; /* f is negative gradient, thus the sign */
}
}
/* Sum the gradient along the line across CPUs */
*
* If we already found a lower value we just skip this step and continue to the update.
*/
+ double gpb;
if (!foundlower)
{
nminstep = 0;
}
/* Take a trial step to this new point - new coords in s_b */
- do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, b, s_min->s.cg_p, s_b,
+ do_em_step(cr, inputrec, mdatoms, fr->bMolPBC, s_min, b, &s_min->s.cg_p, s_b,
constr, top, nrnb, wcycle, -1);
neval++;
/* p does not change within a step, but since the domain decomposition
* might change, we have to use cg_p of s_b here.
*/
- p = s_b->s.cg_p;
- sf = s_b->f;
- gpb = 0;
- for (i = 0; i < mdatoms->homenr; i++)
+ const rvec *pb = as_rvec_array(s_b->s.cg_p.data());
+ const rvec *sfb = as_rvec_array(s_b->f.data());
+ gpb = 0;
+ for (int i = 0; i < mdatoms->homenr; i++)
{
for (m = 0; m < DIM; m++)
{
- gpb -= p[i][m]*sf[i][m]; /* f is negative gradient, thus the sign */
+ gpb -= pb[i][m]*sfb[i][m]; /* f is negative gradient, thus the sign */
}
}
/* Sum the gradient along the line across CPUs */
if (gpb > 0)
{
/* Replace c endpoint with b */
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
c = b;
gpc = gpb;
}
else
{
/* Replace a endpoint with b */
- swap_em_state(s_b, s_a);
+ swap_em_state(&s_b, &s_a);
a = b;
gpa = gpb;
}
fprintf(debug, "CGE: C (%f) is lower than A (%f), moving C to B\n",
s_c->epot, s_a->epot);
}
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
gpb = gpc;
}
else
fprintf(debug, "CGE: A (%f) is lower than C (%f), moving A to B\n",
s_a->epot, s_c->epot);
}
- swap_em_state(s_b, s_a);
+ swap_em_state(&s_b, &s_a);
gpb = gpa;
}
fprintf(debug, "CGE: Found a lower energy %f, moving C to B\n",
s_c->epot);
}
- swap_em_state(s_b, s_c);
+ swap_em_state(&s_b, &s_c);
gpb = gpc;
}
/* update positions */
- swap_em_state(s_min, s_b);
+ swap_em_state(&s_min, &s_b);
gpa = gpb;
/* Print it if necessary */
/* Store the new (lower) energies */
upd_mdebin(mdebin, FALSE, FALSE, (double)step,
mdatoms->tmass, enerd, &s_min->s, inputrec->fepvals, inputrec->expandedvals, s_min->s.box,
- NULL, NULL, vir, pres, NULL, mu_tot, constr);
+ nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
do_log = do_per_step(step, inputrec->nstlog);
do_ene = do_per_step(step, inputrec->nstenergy);
print_ebin_header(fplog, step, step);
}
print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
- do_log ? fplog : NULL, step, step, eprNORMAL,
+ do_log ? fplog : nullptr, step, step, eprNORMAL,
mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
/* Send energies and positions to the IMD client if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
{
/* Write final energy file entries */
print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
- !do_log ? fplog : NULL, step, step, eprNORMAL,
+ !do_log ? fplog : nullptr, step, step, eprNORMAL,
mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
}
write_em_traj(fplog, cr, outf, do_x, do_f, ftp2fn(efSTO, nfile, fnm),
top_global, inputrec, step,
- s_min, state_global);
+ s_min, state_global, energyHistory);
if (MASTER(cr))
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
- fnormn = s_min->fnorm/sqrtNumAtoms;
print_converged(stderr, CG, inputrec->em_tol, step, converged, number_steps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
print_converged(fplog, CG, inputrec->em_tol, step, converged, number_steps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
}
- finish_em(cr, outf, walltime_accounting, wcycle);
+ finish_em(cr, outf, inputrec, walltime_accounting, wcycle);
/* To print the actual number of steps we needed somewhere */
walltime_accounting_set_nsteps_done(walltime_accounting, step);
/*! \brief Do L-BFGS conjugate gradients minimization
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
- int nfile, const t_filenm fnm[],
- const gmx_output_env_t *oenv, gmx_bool bVerbose,
- int nstglobalcomm,
- gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout,
- t_inputrec *inputrec,
- gmx_mtop_t *top_global, t_fcdata *fcd,
- t_state *state_global,
- t_mdatoms *mdatoms,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed,
- t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
- real cpt_period, real max_hours,
- int imdport,
- unsigned long Flags,
- gmx_walltime_accounting_t walltime_accounting)
+ \copydoc integrator_t(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+ int nfile, const t_filenm fnm[],
+ const gmx_output_env_t *oenv, gmx_bool bVerbose,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout,
+ t_inputrec *inputrec,
+ gmx_mtop_t *top_global, t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed,
+ t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ real cpt_period, real max_hours,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
*/
-double do_lbfgs(FILE *fplog, t_commrec *cr,
+double do_lbfgs(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t gmx_unused *oenv, gmx_bool bVerbose,
int gmx_unused nstglobalcomm,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
+ energyhistory_t *energyHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_edsam_t gmx_unused ed,
em_state_t ems;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
int ncorr, nmaxcorr, point, cp, neval, nminstep;
double stepsize, step_taken, gpa, gpb, gpc, tmp, minstep;
- real *rho, *alpha, *ff, *xx, *p, *s, *lastx, *lastf, **dx, **dg;
- real *xa, *xb, *xc, *fa, *fb, *fc, *xtmp, *ftmp;
+ real *rho, *alpha, *p, *s, **dx, **dg;
real a, b, c, maxdelta, delta;
- real diag, Epot0, Epot, EpotA, EpotB, EpotC;
+ real diag, Epot0;
real dgdx, dgdg, sq, yr, beta;
t_mdebin *mdebin;
gmx_bool converged;
rvec mu_tot;
- real fnorm, fmax;
gmx_bool do_log, do_ene, do_x, do_f, foundlower, *frozen;
tensor vir, pres;
int start, end, number_steps;
gmx_mdoutf_t outf;
- int i, k, m, n, nfmax, gf, step;
+ int i, k, m, n, gf, step;
int mdof_flags;
if (PAR(cr))
gmx_fatal(FARGS, "Cannot do parallel L-BFGS Minimization - yet.\n");
}
- if (NULL != constr)
+ if (nullptr != constr)
{
gmx_fatal(FARGS, "The combination of constraints and L-BFGS minimization is not implemented. Either do not use constraints, or use another minimizer (e.g. steepest descent).");
}
n = 3*state_global->natoms;
nmaxcorr = inputrec->nbfgscorr;
- /* Allocate memory */
- /* Use pointers to real so we dont have to loop over both atoms and
- * dimensions all the time...
- * x/f are allocated as rvec *, so make new x0/f0 pointers-to-real
- * that point to the same memory.
- */
- snew(xa, n);
- snew(xb, n);
- snew(xc, n);
- snew(fa, n);
- snew(fb, n);
- snew(fc, n);
snew(frozen, n);
snew(p, n);
- snew(lastx, n);
- snew(lastf, n);
snew(rho, nmaxcorr);
snew(alpha, nmaxcorr);
/* Init em */
init_em(fplog, LBFGS, cr, inputrec,
- state_global, top_global, &ems, &top, &f,
- nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
+ state_global, top_global, &ems, &top,
+ nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
+ vsite, constr, nullptr,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
- /* Do_lbfgs is not completely updated like do_steep and do_cg,
- * so we free some memory again.
- */
- sfree(ems.s.x);
- sfree(ems.f);
-
- xx = (real *)state_global->x;
- ff = (real *)f;
start = 0;
end = mdatoms->homenr;
+ /* We need 4 working states */
+ em_state_t s0 {}, s1 {}, s2 {}, s3 {};
+ em_state_t *sa = &s0;
+ em_state_t *sb = &s1;
+ em_state_t *sc = &s2;
+ em_state_t *last = &s3;
+ /* Initialize by copying the state from ems (we could skip x and f here) */
+ *sa = ems;
+ *sb = ems;
+ *sc = ems;
+
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, LBFGS);
if (vsite)
{
- construct_vsites(vsite, state_global->x, 1, NULL,
+ construct_vsites(vsite, as_rvec_array(state_global->x.data()), 1, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state_global->box);
}
* We do not unshift, so molecules are always whole
*/
neval++;
- ems.s.x = state_global->x;
- ems.f = f;
evaluate_energy(fplog, cr,
top_global, &ems, top,
inputrec, nrnb, wcycle, gstat,
/* Copy stuff to the energy bin for easy printing etc. */
upd_mdebin(mdebin, FALSE, FALSE, (double)step,
mdatoms->tmass, enerd, state_global, inputrec->fepvals, inputrec->expandedvals, state_global->box,
- NULL, NULL, vir, pres, NULL, mu_tot, constr);
+ nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
print_ebin_header(fplog, step, step);
print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
}
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
* vector (force vector the first time), we scale it by the
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
fprintf(stderr, "Using %d BFGS correction steps.\n\n", nmaxcorr);
- fprintf(stderr, " F-max = %12.5e on atom %d\n", fmax, nfmax+1);
- fprintf(stderr, " F-Norm = %12.5e\n", fnorm/sqrtNumAtoms);
+ fprintf(stderr, " F-max = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
+ fprintf(stderr, " F-Norm = %12.5e\n", ems.fnorm/sqrtNumAtoms);
fprintf(stderr, "\n");
/* and copy to the log file too... */
fprintf(fplog, "Using %d BFGS correction steps.\n\n", nmaxcorr);
- fprintf(fplog, " F-max = %12.5e on atom %d\n", fmax, nfmax+1);
- fprintf(fplog, " F-Norm = %12.5e\n", fnorm/sqrtNumAtoms);
+ fprintf(fplog, " F-max = %12.5e on atom %d\n", ems.fmax, ems.a_fmax + 1);
+ fprintf(fplog, " F-Norm = %12.5e\n", ems.fnorm/sqrtNumAtoms);
fprintf(fplog, "\n");
}
point = 0;
// Set initial search direction to the force (-gradient), or 0 for frozen particles.
+ real *fInit = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
for (i = 0; i < n; i++)
{
if (!frozen[i])
{
- dx[point][i] = ff[i]; /* Initial search direction */
+ dx[point][i] = fInit[i]; /* Initial search direction */
}
else
{
// (the main efficiency in the algorithm comes from changing directions), but
// we still need an initial value, so estimate it as the inverse of the norm
// so we take small steps where the potential fluctuates a lot.
- stepsize = 1.0/fnorm;
+ stepsize = 1.0/ems.fnorm;
/* Start the loop over BFGS steps.
* Each successful step is counted, and we continue until
}
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
- top_global, step, (real)step, state_global, state_global, f);
+ top_global, step, (real)step, &ems.s, state_global, energyHistory, &ems.f);
/* Do the linesearching in the direction dx[point][0..(n-1)] */
/* make s a pointer to current search direction - point=0 first time we get here */
s = dx[point];
+ real *xx = static_cast<real *>(as_rvec_array(ems.s.x.data())[0]);
+ real *ff = static_cast<real *>(as_rvec_array(ems.f.data())[0]);
+
// calculate line gradient in position A
for (gpa = 0, i = 0; i < n; i++)
{
}
// Before taking any steps along the line, store the old position
- for (i = 0; i < n; i++)
- {
- lastx[i] = xx[i];
- lastf[i] = ff[i];
- }
- Epot0 = Epot;
+ *last = ems;
+ real *lastx = static_cast<real *>(as_rvec_array(last->s.x.data())[0]);
+ real *lastf = static_cast<real *>(as_rvec_array(last->f.data())[0]);
+ Epot0 = ems.epot;
- for (i = 0; i < n; i++)
- {
- xa[i] = xx[i];
- }
+ *sa = ems;
/* Take a step downhill.
* In theory, we should find the actual minimum of the function in this
// State "A" is the first position along the line.
// reference position along line is initially zero
- EpotA = Epot0;
a = 0.0;
// Check stepsize first. We do not allow displacements
while (maxdelta > inputrec->em_stepsize);
// Take a trial step and move the coordinate array xc[] to position C
+ real *xc = static_cast<real *>(as_rvec_array(sc->s.x.data())[0]);
for (i = 0; i < n; i++)
{
xc[i] = lastx[i] + c*s[i];
neval++;
// Calculate energy for the trial step in position C
- ems.s.x = (rvec *)xc;
- ems.f = (rvec *)fc;
evaluate_energy(fplog, cr,
- top_global, &ems, top,
+ top_global, sc, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, step, FALSE);
- EpotC = ems.epot;
// Calc line gradient in position C
+ real *fc = static_cast<real *>(as_rvec_array(sc->f.data())[0]);
for (gpc = 0, i = 0; i < n; i++)
{
gpc -= s[i]*fc[i]; /* f is negative gradient, thus the sign */
// This is the max amount of increase in energy we tolerate.
// By allowing VERY small changes (close to numerical precision) we
// frequently find even better (lower) final energies.
- tmp = sqrt(GMX_REAL_EPS)*fabs(EpotA);
+ tmp = sqrt(GMX_REAL_EPS)*fabs(sa->epot);
// Accept the step if the energy is lower in the new position C (compared to A),
// or if it is not significantly higher and the line derivative is still negative.
- if (EpotC < EpotA || (gpc < 0 && EpotC < (EpotA+tmp)))
+ if (sc->epot < sa->epot || (gpc < 0 && sc->epot < (sa->epot + tmp)))
{
// Great, we found a better energy. We no longer try to alter the
// stepsize, but simply accept this new better position. The we select a new
// to take smaller steps along a line. Set fnorm based on the new C position,
// which will be used to update the stepsize to 1/fnorm further down.
foundlower = TRUE;
- fnorm = ems.fnorm;
}
else
{
// I also have a safeguard for potentially really pathological functions so we never
// take more than 20 steps before we give up.
// If we already found a lower value we just skip this step and continue to the update.
- nminstep = 0;
+ real fnorm = 0;
+ nminstep = 0;
do
{
// Select a new trial point B in the interval [A,C].
}
// Take a trial step to point B
+ real *xb = static_cast<real *>(as_rvec_array(sb->s.x.data())[0]);
for (i = 0; i < n; i++)
{
xb[i] = lastx[i] + b*s[i];
neval++;
// Calculate energy for the trial step in point B
- ems.s.x = (rvec *)xb;
- ems.f = (rvec *)fb;
evaluate_energy(fplog, cr,
- top_global, &ems, top,
+ top_global, sb, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, step, FALSE);
- EpotB = ems.epot;
- fnorm = ems.fnorm;
+ fnorm = sb->fnorm;
// Calculate gradient in point B
+ real *fb = static_cast<real *>(as_rvec_array(sb->f.data())[0]);
for (gpb = 0, i = 0; i < n; i++)
{
gpb -= s[i]*fb[i]; /* f is negative gradient, thus the sign */
if (gpb > 0)
{
/* Replace c endpoint with b */
- EpotC = EpotB;
- c = b;
- gpc = gpb;
- /* swap coord pointers b/c */
- xtmp = xb;
- ftmp = fb;
- xb = xc;
- fb = fc;
- xc = xtmp;
- fc = ftmp;
+ c = b;
+ /* swap states b and c */
+ swap_em_state(&sb, &sc);
}
else
{
/* Replace a endpoint with b */
- EpotA = EpotB;
- a = b;
- gpa = gpb;
- /* swap coord pointers a/b */
- xtmp = xb;
- ftmp = fb;
- xb = xa;
- fb = fa;
- xa = xtmp;
- fa = ftmp;
+ a = b;
+ /* swap states a and b */
+ swap_em_state(&sa, &sb);
}
/*
*/
nminstep++;
}
- while ((EpotB > EpotA || EpotB > EpotC) && (nminstep < 20));
+ while ((sb->epot > sa->epot || sb->epot > sc->epot) && (nminstep < 20));
- if (fabs(EpotB-Epot0) < GMX_REAL_EPS || nminstep >= 20)
+ if (fabs(sb->epot - Epot0) < GMX_REAL_EPS || nminstep >= 20)
{
/* OK. We couldn't find a significantly lower energy.
* If ncorr==0 this was steepest descent, and then we give up.
/* Select min energy state of A & C, put the best in xx/ff/Epot
*/
- if (EpotC < EpotA)
+ if (sc->epot < sa->epot)
{
- Epot = EpotC;
/* Use state C */
- for (i = 0; i < n; i++)
- {
- xx[i] = xc[i];
- ff[i] = fc[i];
- }
+ ems = *sc;
step_taken = c;
}
else
{
- Epot = EpotA;
/* Use state A */
- for (i = 0; i < n; i++)
- {
- xx[i] = xa[i];
- ff[i] = fa[i];
- }
+ ems = *sa;
step_taken = a;
}
else
{
/* found lower */
- Epot = EpotC;
/* Use state C */
- for (i = 0; i < n; i++)
- {
- xx[i] = xc[i];
- ff[i] = fc[i];
- }
+ ems = *sc;
step_taken = c;
}
}
}
- /* 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))
{
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
fprintf(stderr, "\rStep %d, Epot=%12.6e, Fnorm=%9.3e, Fmax=%9.3e (atom %d)\n",
- step, Epot, fnorm/sqrtNumAtoms, fmax, nfmax+1);
+ step, ems.epot, ems.fnorm/sqrtNumAtoms, ems.fmax, ems.a_fmax + 1);
fflush(stderr);
}
/* Store the new (lower) energies */
upd_mdebin(mdebin, FALSE, FALSE, (double)step,
mdatoms->tmass, enerd, state_global, inputrec->fepvals, inputrec->expandedvals, state_global->box,
- NULL, NULL, vir, pres, NULL, mu_tot, constr);
+ nullptr, nullptr, vir, pres, nullptr, 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);
}
print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
- do_log ? fplog : NULL, step, step, eprNORMAL,
+ do_log ? fplog : nullptr, step, step, eprNORMAL,
mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
/* Send x and E to IMD client, if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
// Reset stepsize in we are doing more iterations
- stepsize = 1.0/fnorm;
+ stepsize = 1.0/ems.fnorm;
/* 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);
+ converged = converged || (ems.fmax < inputrec->em_tol);
} /* End of the loop */
step--; /* we never took that last step in this case */
}
- if (fmax > inputrec->em_tol)
+ if (ems.fmax > inputrec->em_tol)
{
if (MASTER(cr))
{
if (!do_ene || !do_log) /* Write final energy file entries */
{
print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
- !do_log ? fplog : NULL, step, step, eprNORMAL,
+ !do_log ? fplog : nullptr, step, step, eprNORMAL,
mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
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_global);
+ &ems, state_global, energyHistory);
if (MASTER(cr))
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
print_converged(stderr, LBFGS, inputrec->em_tol, step, converged,
- number_steps, Epot, fmax, nfmax, fnorm/sqrtNumAtoms);
+ number_steps, &ems, sqrtNumAtoms);
print_converged(fplog, LBFGS, inputrec->em_tol, step, converged,
- number_steps, Epot, fmax, nfmax, fnorm/sqrtNumAtoms);
+ number_steps, &ems, sqrtNumAtoms);
fprintf(fplog, "\nPerformed %d energy evaluations in total.\n", neval);
}
- finish_em(cr, outf, walltime_accounting, wcycle);
+ finish_em(cr, outf, inputrec, walltime_accounting, wcycle);
/* To print the actual number of steps we needed somewhere */
walltime_accounting_set_nsteps_done(walltime_accounting, step);
} /* That's all folks */
/*! \brief Do steepest descents minimization
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
- int nfile, const t_filenm fnm[],
- const gmx_output_env_t *oenv, gmx_bool bVerbose,
- int nstglobalcomm,
- gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout,
- t_inputrec *inputrec,
- gmx_mtop_t *top_global, t_fcdata *fcd,
- t_state *state_global,
- t_mdatoms *mdatoms,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed,
- t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
- real cpt_period, real max_hours,
- int imdport,
- unsigned long Flags,
- gmx_walltime_accounting_t walltime_accounting)
+ \copydoc integrator_t(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+ int nfile, const t_filenm fnm[],
+ const gmx_output_env_t *oenv, gmx_bool bVerbose,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout,
+ t_inputrec *inputrec,
+ gmx_mtop_t *top_global, t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed,
+ t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ real cpt_period, real max_hours,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
*/
-double do_steep(FILE *fplog, t_commrec *cr,
+double do_steep(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t gmx_unused *oenv, gmx_bool bVerbose,
int gmx_unused nstglobalcomm,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
+ energyhistory_t *energyHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_edsam_t gmx_unused ed,
gmx_walltime_accounting_t walltime_accounting)
{
const char *SD = "Steepest Descents";
- em_state_t *s_min, *s_try;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
real stepsize;
- real ustep, fnormn;
+ real ustep;
gmx_mdoutf_t outf;
t_mdebin *mdebin;
gmx_bool bDone, bAbort, do_x, do_f;
int count = 0;
int steps_accepted = 0;
- s_min = init_em_state();
- s_try = init_em_state();
+ /* Create 2 states on the stack and extract pointers that we will swap */
+ em_state_t s0 {}, s1 {};
+ em_state_t *s_min = &s0;
+ em_state_t *s_try = &s1;
/* Init em and store the local state in s_try */
init_em(fplog, SD, cr, inputrec,
- state_global, top_global, s_try, &top, &f,
- nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
+ state_global, top_global, s_try, &top,
+ nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
+ vsite, constr, nullptr,
nfile, fnm, &outf, &mdebin, imdport, Flags, wcycle);
/* Print to log file */
{
validStep =
do_em_step(cr, inputrec, mdatoms, fr->bMolPBC,
- s_min, stepsize, s_min->f, s_try,
+ s_min, stepsize, &s_min->f, s_try,
constr, top, nrnb, wcycle, count);
}
/* Store the new (lower) energies */
upd_mdebin(mdebin, FALSE, FALSE, (double)count,
mdatoms->tmass, enerd, &s_try->s, inputrec->fepvals, inputrec->expandedvals,
- s_try->s.box, NULL, NULL, vir, pres, NULL, mu_tot, constr);
+ s_try->s.box, nullptr, nullptr, vir, pres, nullptr, mu_tot, constr);
/* Prepare IMD energy record, if bIMD is TRUE. */
IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, count, TRUE);
/* Copy the arrays for force, positions and energy */
/* The 'Min' array always holds the coords and forces of the minimal
sampled energy */
- swap_em_state(s_min, s_try);
+ swap_em_state(&s_min, &s_try);
if (count > 0)
{
ustep *= 1.2;
/* Write to trn, if necessary */
do_x = do_per_step(steps_accepted, inputrec->nstxout);
do_f = do_per_step(steps_accepted, inputrec->nstfout);
- write_em_traj(fplog, cr, outf, do_x, do_f, NULL,
+ write_em_traj(fplog, cr, outf, do_x, do_f, nullptr,
top_global, inputrec, count,
- s_min, state_global);
+ s_min, state_global, energyHistory);
}
else
{
{
if (MASTER(cr))
{
- warn_step(stderr, inputrec->em_tol, count == nsteps, constr != NULL);
- warn_step(fplog, inputrec->em_tol, count == nsteps, constr != NULL);
+ warn_step(stderr, inputrec->em_tol, count == nsteps, constr != nullptr);
+ warn_step(fplog, inputrec->em_tol, count == nsteps, constr != nullptr);
}
bAbort = TRUE;
}
/* Send IMD energies and positions, if bIMD is TRUE. */
- if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, as_rvec_array(state_global->x.data()), inputrec, 0, wcycle) && MASTER(cr))
{
IMD_send_positions(inputrec->imd);
}
}
write_em_traj(fplog, cr, outf, TRUE, inputrec->nstfout, ftp2fn(efSTO, nfile, fnm),
top_global, inputrec, count,
- s_min, state_global);
+ s_min, state_global, energyHistory);
if (MASTER(cr))
{
double sqrtNumAtoms = sqrt(static_cast<double>(state_global->natoms));
- fnormn = s_min->fnorm/sqrtNumAtoms;
print_converged(stderr, SD, inputrec->em_tol, count, bDone, nsteps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
print_converged(fplog, SD, inputrec->em_tol, count, bDone, nsteps,
- s_min->epot, s_min->fmax, s_min->a_fmax, fnormn);
+ s_min, sqrtNumAtoms);
}
- finish_em(cr, outf, walltime_accounting, wcycle);
+ finish_em(cr, outf, inputrec, walltime_accounting, wcycle);
/* To print the actual number of steps we needed somewhere */
inputrec->nsteps = count;
} /* That's all folks */
/*! \brief Do normal modes analysis
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
- int nfile, const t_filenm fnm[],
- const gmx_output_env_t *oenv, gmx_bool bVerbose,
- int nstglobalcomm,
- gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout,
- t_inputrec *inputrec,
- gmx_mtop_t *top_global, t_fcdata *fcd,
- t_state *state_global,
- t_mdatoms *mdatoms,
- t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed,
- t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
- real cpt_period, real max_hours,
- int imdport,
- unsigned long Flags,
- gmx_walltime_accounting_t walltime_accounting)
+ \copydoc integrator_t(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+ int nfile, const t_filenm fnm[],
+ const gmx_output_env_t *oenv, gmx_bool bVerbose,
+ int nstglobalcomm,
+ gmx_vsite_t *vsite, gmx_constr_t constr,
+ int stepout,
+ t_inputrec *inputrec,
+ gmx_mtop_t *top_global, t_fcdata *fcd,
+ t_state *state_global,
+ t_mdatoms *mdatoms,
+ t_nrnb *nrnb, gmx_wallcycle_t wcycle,
+ gmx_edsam_t ed,
+ t_forcerec *fr,
+ int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ real cpt_period, real max_hours,
+ int imdport,
+ unsigned long Flags,
+ gmx_walltime_accounting_t walltime_accounting)
*/
-double do_nm(FILE *fplog, t_commrec *cr,
+double do_nm(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t gmx_unused *oenv, gmx_bool bVerbose,
int gmx_unused nstglobalcomm,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
+ energyhistory_t gmx_unused *energyHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
gmx_edsam_t gmx_unused ed,
int nnodes, node;
gmx_localtop_t *top;
gmx_enerdata_t *enerd;
- rvec *f;
gmx_global_stat_t gstat;
t_graph *graph;
tensor vir, pres;
rvec *fneg, *dfdx;
gmx_bool bSparse; /* use sparse matrix storage format */
size_t sz;
- gmx_sparsematrix_t * sparse_matrix = NULL;
- real * full_matrix = NULL;
- em_state_t * state_work;
+ gmx_sparsematrix_t * sparse_matrix = nullptr;
+ real * full_matrix = nullptr;
/* added with respect to mdrun */
int row, col;
real x_min;
bool bIsMaster = MASTER(cr);
- if (constr != NULL)
+ if (constr != nullptr)
{
gmx_fatal(FARGS, "Constraints present with Normal Mode Analysis, this combination is not supported");
}
- state_work = init_em_state();
+ gmx_shellfc_t *shellfc;
+
+ em_state_t state_work {};
/* Init em and store the local state in state_minimum */
init_em(fplog, NM, cr, inputrec,
- state_global, top_global, state_work, &top,
- &f,
- nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
- nfile, fnm, &outf, NULL, imdport, Flags, wcycle);
-
- gmx_shellfc_t *shellfc = init_shell_flexcon(stdout,
- top_global,
- n_flexible_constraints(constr),
- inputrec->nstcalcenergy,
- DOMAINDECOMP(cr));
+ state_global, top_global, &state_work, &top,
+ nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat,
+ vsite, constr, &shellfc,
+ nfile, fnm, &outf, nullptr, imdport, Flags, wcycle);
- if (shellfc)
- {
- make_local_shells(cr, mdatoms, shellfc);
- }
std::vector<size_t> atom_index = get_atom_index(top_global);
snew(fneg, atom_index.size());
snew(dfdx, atom_index.size());
*/
if (EEL_FULL(fr->eeltype) || fr->rlist == 0.0)
{
- md_print_info(cr, fplog, "Non-cutoff electrostatics used, forcing full Hessian format.\n");
+ GMX_LOG(mdlog.warning).appendText("Non-cutoff electrostatics used, forcing full Hessian format.");
bSparse = FALSE;
}
else if (atom_index.size() < 1000)
{
- md_print_info(cr, fplog, "Small system size (N=%d), using full Hessian format.\n", atom_index.size());
+ GMX_LOG(mdlog.warning).appendTextFormatted("Small system size (N=%d), using full Hessian format.",
+ atom_index.size());
bSparse = FALSE;
}
else
{
- md_print_info(cr, fplog, "Using compressed symmetric sparse Hessian format.\n");
+ GMX_LOG(mdlog.warning).appendText("Using compressed symmetric sparse Hessian format.");
bSparse = TRUE;
}
/* Make evaluate_energy do a single node force calculation */
cr->nnodes = 1;
evaluate_energy(fplog, cr,
- top_global, state_work, top,
+ top_global, &state_work, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, -1, TRUE);
cr->nnodes = nnodes;
/* if forces are not small, warn user */
- get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, state_work);
+ get_state_f_norm_max(cr, &(inputrec->opts), mdatoms, &state_work);
- md_print_info(cr, fplog, "Maximum force:%12.5e\n", state_work->fmax);
- if (state_work->fmax > 1.0e-3)
+ GMX_LOG(mdlog.warning).appendTextFormatted("Maximum force:%12.5e", state_work.fmax);
+ if (state_work.fmax > 1.0e-3)
{
- md_print_info(cr, fplog,
- "The force is probably not small enough to "
- "ensure that you are at a minimum.\n"
- "Be aware that negative eigenvalues may occur\n"
- "when the resulting matrix is diagonalized.\n\n");
+ GMX_LOG(mdlog.warning).appendText(
+ "The force is probably not small enough to "
+ "ensure that you are at a minimum.\n"
+ "Be aware that negative eigenvalues may occur\n"
+ "when the resulting matrix is diagonalized.");
}
/***********************************************************
int force_flags = GMX_FORCE_STATECHANGED | GMX_FORCE_ALLFORCES;
double t = 0;
- x_min = state_work->s.x[atom][d];
+ x_min = state_work.s.x[atom][d];
for (unsigned int dx = 0; (dx < 2); dx++)
{
if (dx == 0)
{
- state_work->s.x[atom][d] = x_min - der_range;
+ state_work.s.x[atom][d] = x_min - der_range;
}
else
{
- state_work->s.x[atom][d] = x_min + der_range;
+ state_work.s.x[atom][d] = x_min + der_range;
}
/* Make evaluate_energy do a single node force calculation */
inputrec, bNS, force_flags,
top,
constr, enerd, fcd,
- &state_work->s, state_work->f, vir, mdatoms,
+ &state_work.s, &state_work.f, vir, mdatoms,
nrnb, wcycle, graph, &top_global->groups,
shellfc, fr, bBornRadii, t, mu_tot,
- vsite, NULL);
+ vsite);
bNS = false;
step++;
}
else
{
evaluate_energy(fplog, cr,
- top_global, state_work, top,
+ top_global, &state_work, top,
inputrec, nrnb, wcycle, gstat,
vsite, constr, fcd, graph, mdatoms, fr,
mu_tot, enerd, vir, pres, atom*2+dx, FALSE);
{
for (size_t i = 0; i < atom_index.size(); i++)
{
- copy_rvec(state_work->f[atom_index[i]], fneg[i]);
+ copy_rvec(state_work.f[atom_index[i]], fneg[i]);
}
}
}
/* x is restored to original */
- state_work->s.x[atom][d] = x_min;
+ state_work.s.x[atom][d] = x_min;
for (size_t j = 0; j < atom_index.size(); j++)
{
for (size_t k = 0; (k < DIM); k++)
{
dfdx[j][k] =
- -(state_work->f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
+ -(state_work.f[atom_index[j]][k] - fneg[j][k])/(2*der_range);
}
}
gmx_mtxio_write(ftp2fn(efMTX, nfile, fnm), sz, sz, full_matrix, sparse_matrix);
}
- finish_em(cr, outf, walltime_accounting, wcycle);
+ finish_em(cr, outf, inputrec, walltime_accounting, wcycle);
walltime_accounting_set_nsteps_done(walltime_accounting, atom_index.size()*2);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff