diff --git a/.gitignore b/.gitignore index 95327c32..631e65ab 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,8 @@ *.o *.mod *.tgz +*.i90 +*__genmod.f90 github_bin/ build_majestix build_commands diff --git a/CMakeLists.txt b/CMakeLists.txt index 37fc48df..d8da2cae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,131 +10,260 @@ # 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 crest. If not, see . + +cmake_minimum_required(VERSION 3.17) +option(WITH_OBJECT "To build using object library" TRUE) +option(INSTALL_MODULES "Install Fortran module files to include directory." FALSE) -cmake_minimum_required(VERSION 3.14) +# Buggy CMake versions +if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.27.0 AND CMAKE_VERSION VERSION_LESS 3.28.0) + set(WITH_OBJECT FALSE) +endif() +# Setup the crest Project project( - "crest" - LANGUAGES "Fortran" "C" - VERSION "3.0" + crest + LANGUAGES "C" "Fortran" + VERSION 3.0 + DESCRIPTION "A tool for the exploration of low-energy chemical space" ) set(SOVERSION "pre") +#enable_testing() + # Follow GNU conventions for installing directories include(GNUInstallDirs) -# General configuration information -set(libs) +# Include further configurations add_subdirectory("config") -# -# Libraries -# +############################################################################### +####################### SUBPROJECTS & DEPENDENCIES ############################ +############################################################################### +# Check a specific CMake targets and execute the +# corresponding Find scripts (located in config/modules/) -# OpenMP dependency -if(NOT TARGET OpenMP::OpenMP_Fortran AND WITH_OpenMP) - find_package(OpenMP REQUIRED) - list( - APPEND libs - OpenMP::OpenMP_Fortran - ) +# LAPACK, BLAS, OpenMP (usually via shared dependencies) +if(STATICBUILD) + set(BLA_STATIC ON) endif() - -# BLAS and LAPACK -if(NOT TARGET BLAS::BLAS) - find_package(BLAS REQUIRED) - if(NOT TARGET BLAS::BLAS AND BLAS_FOUND) - add_library(BLAS::BLAS INTERFACE IMPORTED) - target_link_libraries(BLAS::BLAS INTERFACE "${BLAS_LIBRARIES}") - target_link_options(BLAS::BLAS INTERFACE "${BLAS_LINKER_FLAGS}") - endif() +find_package(LAPACK REQUIRED) +find_package(BLAS REQUIRED) +if(NOT TARGET "OpenMP::OpenMP_Fortran" AND WITH_OpenMP) + find_package("OpenMP" REQUIRED) endif() -if(NOT TARGET LAPACK::LAPACK) - find_package(LAPACK REQUIRED) - if(NOT TARGET LAPACK::LAPACK AND LAPACK_FOUND) - add_library(LAPACK::LAPACK INTERFACE IMPORTED) - target_link_libraries(LAPACK::LAPACK INTERFACE "${LAPACK_LIBRARIES}") - target_link_options(LAPACK::LAPACK INTERFACE "${LAPACK_LINKER_FLAGS}") - endif() + +# TOML-F (needs to be imported before tblite) +if(NOT TARGET "toml-f::toml-f" AND WITH_TOMLF) + find_package("toml-f" REQUIRED) + add_compile_definitions(WITH_TOMLF) endif() -list( - APPEND libs - LAPACK::LAPACK - BLAS::BLAS -) # tblite -if(WITH_TBLITE) - find_package("tblite" REQUIRED) - add_definitions(-DWITH_TBLITE) - list( - APPEND libs - tblite::tblite - ) +if(NOT TARGET "tblite::tblite" AND WITH_TBLITE) + find_package("tblite" REQUIRED) + add_compile_definitions(WITH_TBLITE) endif() +# GFN-FF +if(NOT TARGET "gfnff::gfnff" AND WITH_GFNFF) + find_package("gfnff" REQUIRED) + add_compile_definitions(WITH_GFNFF) +endif() -# toml-f -if(WITH_TOMLF) - find_package("toml-f" REQUIRED) - add_definitions(-DWITH_TOMLF) - list( - APPEND libs - toml-f::toml-f - ) +# GFN0-xTB +if(NOT TARGET "gfn0::gfn0" AND WITH_GFN0) + find_package("gfn0" REQUIRED) + add_compile_definitions(WITH_GFN0) endif() +# XHCFF +if(NOT TARGET "xhcff::xhcff" AND WITH_XHCFF) + find_package("xhcff" REQUIRED) + add_compile_definitions(WITH_XHCFF) +endif() -# GFN0-xTB -if(WITH_GFN0) - find_package("gfn0" REQUIRED) - add_definitions(-DWITH_GFN0) - list( - APPEND libs - gfn0::gfn0 - ) +# lwONIOM +if(NOT TARGET "lwoniom::lwoniom" AND WITH_LWONIOM) + find_package("lwoniom" REQUIRED) + add_compile_definitions(WITH_LWONIOM) endif() +# Sources: initialize program sources (prog) and library sources (srcs) empty +set(prog) +set(srcs) +add_subdirectory("src") -# GFN-FF -if(WITH_GFNFF) - find_package("gfnff" REQUIRED) - add_definitions(-DWITH_GFNFF) - list( - APPEND libs - gfnff::gfnff - ) + +if(NOT EXISTS "${PROJECT_BINARY_DIR}/include") + file(MAKE_DIRECTORY "${PROJECT_BINARY_DIR}/include") endif() -# XHCFF -if(WITH_XHCFF) - find_package("xhcff" REQUIRED) - add_definitions(-DWITH_XHCFF) - list( - APPEND libs - xhcff::xhcff - ) +############################################################################### +########################## OBJECT LIBRARY ##################################### +############################################################################### +if(WITH_OBJECT AND NOT STATICBUILD) + + add_library( + "${PROJECT_NAME}-object" + OBJECT + ${srcs} + ) + + # customize object library + set_target_properties( + "${PROJECT_NAME}-object" + PROPERTIES + Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include" + POSITION_INDEPENDENT_CODE ON + ) + + # link object library conditionally against required + target_link_libraries( + "${PROJECT_NAME}-object" + PUBLIC + $<$:tblite::tblite> + $<$:gfn0::gfn0> + $<$:gfnff::gfnff> + $<$:xhcff::xhcff> + $<$:toml-f::toml-f> + $<$:lwoniom::lwoniom> + $<$:OpenMP::OpenMP_Fortran> + ) + + # include directories + target_include_directories( + "${PROJECT_NAME}-object" + PUBLIC + ${crest-config-dir} + ${PROJECT_BINARY_DIR} + $ + $ + $/${CMAKE_INSTALL_INCLUDEDIR}> + ) + endif() +############################################################################### +############################### Static Library ################################ +############################################################################### +if(WITH_OBJECT AND NOT STATICBUILD) + add_library( + "lib-${PROJECT_NAME}-static" + STATIC + $ + ) +else() + add_library( + "lib-${PROJECT_NAME}-static" + STATIC + ${srcs} + ) +endif() +# workaround attemp. doesn't work yet +set(LINK_OpenMP FALSE) +if(WITH_OpenMP AND NOT STATICBUILD) + set(LINK_OpenMP TRUE) +endif() +target_link_libraries( + "lib-${PROJECT_NAME}-static" + PUBLIC + ${BLAS_LIBRARIES} + ${LAPACK_LIBRARIES} + $<$:OpenMP::OpenMP_Fortran> + $<$:tblite::tblite> + $<$:gfn0::gfn0> + $<$:gfnff::gfnff> + $<$:xhcff::xhcff> + $<$:toml-f::toml-f> + $<$:lwoniom::lwoniom> + $<$:-static> +) -# -# CREST sources -# -set(prog) -set(srcs) -add_subdirectory("src") +set_target_properties( + "lib-${PROJECT_NAME}-static" + PROPERTIES + Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include" + ARCHIVE_OUTPUT_DIRECTORY "${PROJECT_BINARY_DIR}" + POSITION_INDEPENDENT_CODE ON + OUTPUT_NAME "${PROJECT_NAME}" + ) -# -# Executables -# + +target_include_directories( + "lib-${PROJECT_NAME}-static" + PUBLIC + $ + $ + $/${CMAKE_INSTALL_INCLUDEDIR}> +) + + +#################################################################################### +############################# Shared Library ####################################### +#################################################################################### +if (WITH_OBJECT AND NOT STATICBUILD) + add_library( + "lib-${PROJECT_NAME}-shared" + SHARED + $ + ) + + target_link_libraries( + "lib-${PROJECT_NAME}-shared" + PUBLIC + ${BLAS_LIBRARIES} + ${LAPACK_LIBRARIES} + $<$:OpenMP::OpenMP_Fortran> + $<$:tblite::tblite> + $<$:gfn0::gfn0> + $<$:gfnff::gfnff> + $<$:xhcff::xhcff> + $<$:toml-f::toml-f> + $<$:lwoniom::lwoniom> + ) + + set_target_properties( + "lib-${PROJECT_NAME}-shared" + PROPERTIES + Fortran_MODULE_DIRECTORY "${PROJECT_BINARY_DIR}/include" + LIBRARY_OUTPUT_DIRECTORY "${PROJECT_BINARY_DIR}" + OUTPUT_NAME "${PROJECT_NAME}" + VERSION "${PROJECT_VERSION}" + SOVERSION "${PROJECT_VERSION_MAJOR}" + ) + + target_include_directories( + "lib-${PROJECT_NAME}-shared" + PUBLIC + $ + $ + $/${CMAKE_INSTALL_INCLUDEDIR}> + ) +endif() + +############################################################################### +############################### Executables ################################### +############################################################################### add_executable( ${PROJECT_NAME}-exe - "${srcs}" "${prog}" + ${prog} ) + +target_link_libraries( + ${PROJECT_NAME}-exe + PRIVATE + "lib-${PROJECT_NAME}-static" + $<$:-static> +) + set_target_properties( ${PROJECT_NAME}-exe PROPERTIES @@ -142,24 +271,11 @@ set_target_properties( RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR} OUTPUT_NAME "${PROJECT_NAME}" ) -target_link_libraries( - ${PROJECT_NAME}-exe - PRIVATE - "${libs}" -) - -get_target_property(OUT ${PROJECT_NAME}-exe LINK_LIBRARIES) -message(STATUS ${OUT}) -# -# Binary installing option -# -install( - TARGETS - "${PROJECT_NAME}-exe" - EXPORT - "${PROJECT_NAME}-targets" - RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" - LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" - ARCHIVE DESTINATION "${CMAKE_INSTALL_LIBDIR}" +target_include_directories( + ${PROJECT_NAME}-exe + PRIVATE + ${PROJECT_SOURCE_DIR}/include + $ ) + diff --git a/README.md b/README.md index 10e009f7..07ad9188 100644 --- a/README.md +++ b/README.md @@ -63,6 +63,8 @@ Working and tested builds of CREST (mostly on Ubuntu 20.04 LTS): |--------------|----------|------------------------|:--------------:|:----------:| | CMake | GNU (gcc 10.3.0) | [OpenBLAS](https://github.com/xianyi/OpenBLAS) (with OpenMP) | dynamic | ✅ | | CMake | GNU (gcc 10.3.0) | [MKL shared (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ✅ | +| CMake | GNU (gcc 9.5.0) | [MKL shared (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ✅ | +| CMake | [Intel (`ifort`/`icc` 2021.9.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html) | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | dynamic | ✅ | | Meson | [Intel (`ifort`/`icc` 2021.9.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html) | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | static | ✅ | | Meson | [Intel (`ifort` 2021.9.0/`icx` 2023.1.0)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html) | [MKL static (oneAPI 2023.1)](https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl.html) | static | ✅ | @@ -77,7 +79,26 @@ For more information about builds including subprojects see [here](./subprojects Some basic build instructions can be found in the following dropdown tabs: + +
+

CMake build

 + + +Building CREST with CMake works with the following chain of commands (in this example with `gfortran/gcc` compilers): +```bash +export FC=gfortran CC=gcc +cmake -B _build -DCMAKE_BUILD_TYPE=Release +``` +and then to build the CREST binary +```bash +make -C _build +``` + +The `CMake` build typically requires access to shared libraries of LAPACK and OpenMP. They must be present in the library paths at compile and runtime. +
+ +

meson build

 @@ -97,25 +118,6 @@ When attempting to build with `gfortran` and `gcc`, add `-Dla_backend=mkl` to th By default the `meson` build will create a **statically** linked binary.
-
-

CMake build

 - - -For the setup of CMake see also the [CMake setup](https://github.com/grimme-lab/xtb/blob/master/cmake/README.adoc) page hosted at the `xtb` repository. -Building CREST with CMake works with the following chain of commands: -```bash -export FC=gfortran CC=gcc -cmake -B _build -DCMAKE_BUILD_TYPE=Release -``` -and then to build the CREST binary -```bash -make -C _build -``` - -The CMake build of CREST is focused on and tested with the GNU `gfortran`/`gcc` compilers. The Intel compilers could technically be used as well, but in our experience the respective build is more fragile than its static `meson` counterpart. - -By default the `CMake` build will create a **dynamically** linked binary. -

Conda build

 diff --git a/assets/template/metadata.f90 b/assets/template/metadata.f90 index b2e4da29..38a4e74b 100644 --- a/assets/template/metadata.f90 +++ b/assets/template/metadata.f90 @@ -10,4 +10,4 @@ character(len=*),parameter :: gfnffvar = "@gfnffvar@" character(len=*),parameter :: tblitevar = "@tblitevar@" character(len=*),parameter :: xhcffvar = "@xhcffvar@" - +character(len=*),parameter :: lwoniomvar = "@lwoniomvar@" diff --git a/config/CMakeLists.txt b/config/CMakeLists.txt index 06497a99..e6a7e09d 100644 --- a/config/CMakeLists.txt +++ b/config/CMakeLists.txt @@ -13,21 +13,34 @@ # # You should have received a copy of the GNU Lesser General Public License # along with crest. If not, see . +# Set the module path for CMake includes +######################################################################################### +######################################################################################### +# Add modules to the CMake build list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/modules") set(CMAKE_MODULE_PATH "${CMAKE_MODULE_PATH}" PARENT_SCOPE) + +# specify module installation directory install( DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/modules/" DESTINATION "${CMAKE_INSTALL_LIBDIR}/cmake/${PROJECT_NAME}" ) -set( - module-dir - "${PROJECT_NAME}/${CMAKE_Fortran_COMPILER_ID}-${CMAKE_Fortran_COMPILER_VERSION}" -) -set(module-dir "${module-dir}" PARENT_SCOPE) +# Options for enabling or disabling features +option(WITH_OpenMP "Enable OpenMP support" TRUE) +option(WITH_TBLITE "Enable support for tblite" TRUE) +option(WITH_TOMLF "Enable support for toml-f" TRUE) +option(WITH_GFN0 "Enable support for GFN0-xTB" TRUE) +option(WITH_GFNFF "Enable support for GFN-FF" TRUE) +option(WITH_XHCFF "Enable support for XHCFF" FALSE) +option(WITH_LWONIOM "Enable support for lwONIOM" TRUE) +option(STATICBUILD "Attempt to link everything statically" FALSE) # doesn't work yet + +######################################################################################### +######################################################################################### # Set build type as CMake does not provide defaults if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) @@ -42,57 +55,28 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) ) endif() +######################################################################################### +######################################################################################### -### Options and defaults - -include("${CMAKE_CURRENT_SOURCE_DIR}/modules/crest-utils.cmake") -set(fortran_minimal_versions "GNU;7.5" "Intel;19.0") -check_minimal_compiler_version("Fortran" "${fortran_minimal_versions}") - -option(WITH_OpenMP "Enable support for shared memory parallelisation with OpenMP" TRUE) - -option(WITH_TBLITE "Enable build with the lightweight tight-binding library" TRUE) - -option(WITH_TOMLF "Enable build with toml-f support" TRUE) - -option(WITH_GFN0 "Enable build with GFN0-xTB support" TRUE) - -option(WITH_GFNFF "Enable build with GFN-FF support" TRUE) - -option(WITH_XHCFF "Enable build with XHCFF support" FALSE) - -option(WITH_LWONIOM "Enable build with lwONIOM support" TRUE) - -if(NOT DEFINED "${PROJECT_NAME}-dependency-method") - set( - "${PROJECT_NAME}-dependency-method" - "subproject" "cmake" "pkgconf" "fetch" - ) -endif() - -# -# Compiler settings -# +# Compiler settings for GNU and Intel Fortran compilers if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") -# set(dialect "-fdefault-real-8 -fdefault-double-8 -ffree-line-length-none -fbacktrace") set(dialect "-g -O0 -fbacktrace -ffree-line-length-none -fbacktrace") - set(bounds "-fbounds-check") -endif() -if(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") -# set(dialect "-axAVX2 -r8 -traceback") - set(dialect "-g -O2 -r8 -align array64byte -traceback") - set(bounds "-check bounds") -endif() -if(CMAKE_Fortran_COMPILER_ID MATCHES "PGI") - set(dialect "-Mbackslash -Mallocatable=03 -r8 -traceback") + set(bounds "-fbounds-check -ffpe-trap=invalid,zero,overflow") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "Intel") + set(dialect "-g -O2 -r8 -align array64byte -traceback") + set(bounds "-check all -fpe0") +else() + message(FATAL_ERROR "Please use an Intel or GNU compiler!") endif() + +# Apply the compiler flags set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} ${bounds}" PARENT_SCOPE) set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} ${dialect}" PARENT_SCOPE) +######################################################################################### +######################################################################################### -# # Populate crest_metadata.fh -# set(version ${PROJECT_VERSION}${SOVERSION}) execute_process(COMMAND git show -s --format=%h RESULT_VARIABLE git_return @@ -111,7 +95,7 @@ set(fcver ${CMAKE_Fortran_COMPILER_VERSION}) set(ccid ${CMAKE_C_COMPILER_ID}) set(ccver ${CMAKE_C_COMPILER_VERSION}) set(bsystem "cmake ${CMAKE_VERSION}") -set(tomlfvar "${WITH_TOMLF}")#string(TOLOWER "${MY_BOOL}" ${WITH_TOMLF})) +set(tomlfvar "${WITH_TOMLF}") set(gfn0var "${WITH_GFN0}") set(gfnffvar "${WITH_GFNFF}") set(tblitevar "${WITH_TBLITE}") @@ -123,12 +107,6 @@ configure_file( "${PROJECT_BINARY_DIR}/crest_metadata.fh" @ONLY ) -# gfortran needs the file at another place... -configure_file( - "${PROJECT_SOURCE_DIR}/assets/template/metadata.f90" - "${PROJECT_BINARY_DIR}/include/crest_metadata.fh" - @ONLY -) - - +######################################################################################### +######################################################################################### diff --git a/config/modules/Findgfn0.cmake b/config/modules/Findgfn0.cmake index f5699b7a..a44f8c8d 100644 --- a/config/modules/Findgfn0.cmake +++ b/config/modules/Findgfn0.cmake @@ -19,29 +19,19 @@ set(_pkg "GFN0") set(_url "https://github.com/pprcht/gfn0") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf") endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") crest_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}") +set (found FALSE) if(TARGET "gfn0::gfn0") set (found TRUE) -else() - set (found FALSE) endif() -message("-- Found GFN0: ${found}") +message(STATUS "Found GFN0-xTB: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/Findgfnff.cmake b/config/modules/Findgfnff.cmake index 9308b613..faf16024 100644 --- a/config/modules/Findgfnff.cmake +++ b/config/modules/Findgfnff.cmake @@ -19,29 +19,19 @@ set(_pkg "GFNFF") set(_url "https://github.com/pprcht/gfnff") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf") endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") crest_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}") +set(found FALSE) if(TARGET "gfnff::gfnff") set (found TRUE) -else() - set (found FALSE) endif() -message("-- Found GFN-FF: ${found}") +message(STATUS "Found GFN-FF: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/Findlwoniom.cmake b/config/modules/Findlwoniom.cmake index 8c131dbe..1f825b13 100644 --- a/config/modules/Findlwoniom.cmake +++ b/config/modules/Findlwoniom.cmake @@ -16,32 +16,22 @@ set(_lib "lwoniom") set(_pkg "LWONIOM") -set(_url "https://github.com/pprcht/lwoniom") +set(_url "https://github.com/crest-lab/lwoniom") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf") endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") crest_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}") +set(found FALSE) if(TARGET "lwoniom::lwoniom") set (found TRUE) -else() - set (found FALSE) endif() -message("-- Found lwONIOM: ${found}") +message(STATUS "Found lwONIOM: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/Findtblite.cmake b/config/modules/Findtblite.cmake index 51464b5d..9ca3cfc3 100644 --- a/config/modules/Findtblite.cmake +++ b/config/modules/Findtblite.cmake @@ -19,29 +19,19 @@ set(_pkg "TBLITE") set(_url "https://github.com/tblite/tblite") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf") endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") crest_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}") +set(found FALSE) if(TARGET "tblite::tblite") set (found TRUE) -else() - set (found FALSE) endif() -message("-- Found tblite: ${found}") +message(STATUS "Found tblite: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/Findtoml-f.cmake b/config/modules/Findtoml-f.cmake index 81eeb159..ca1da4cb 100644 --- a/config/modules/Findtoml-f.cmake +++ b/config/modules/Findtoml-f.cmake @@ -19,12 +19,7 @@ set(_pkg "TOML-F") set(_url "https://github.com/toml-f/toml-f") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf") endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") @@ -36,12 +31,8 @@ if(TARGET "toml-f::toml-f") else() set (found FALSE) endif() -message("-- Found toml-f: ${found}") +message(STATUS "Found toml-f: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/Findxhcff.cmake b/config/modules/Findxhcff.cmake index bcfaae4e..fff680e1 100644 --- a/config/modules/Findxhcff.cmake +++ b/config/modules/Findxhcff.cmake @@ -19,29 +19,19 @@ set(_pkg "XHCFF") set(_url "https://github.com/zellerf/xhcff-lib") if(NOT DEFINED "${_pkg}_FIND_METHOD") - if(DEFINED "${PROJECT_NAME}-dependency-method") - set("${_pkg}_FIND_METHOD" "${${PROJECT_NAME}-dependency-method}") - else() - set("${_pkg}_FIND_METHOD" "cmake" "pkgconf" "subproject" "fetch") - endif() - set("_${_pkg}_FIND_METHOD") + set("${_pkg}_FIND_METHOD" "subproject" "cmake" "fetch" "pkgconf" ) endif() include("${CMAKE_CURRENT_LIST_DIR}/crest-utils.cmake") crest_find_package("${_lib}" "${${_pkg}_FIND_METHOD}" "${_url}") +set(found FALSE) if(TARGET "xhcff::xhcff") set (found TRUE) -else() - set (found FALSE) endif() -message("-- Found xhcff: ${found}") +message(STATUS "Found xhcff: ${found}") -if(DEFINED "_${_pkg}_FIND_METHOD") - unset("${_pkg}_FIND_METHOD") - unset("_${_pkg}_FIND_METHOD") -endif() unset(_lib) unset(_pkg) unset(_url) diff --git a/config/modules/crest-utils.cmake b/config/modules/crest-utils.cmake index dfc38046..375714a2 100644 --- a/config/modules/crest-utils.cmake +++ b/config/modules/crest-utils.cmake @@ -14,6 +14,9 @@ # You should have received a copy of the GNU Lesser General Public License # along with crest. If not, see . +######################################################################################### +######################################################################################### + # Handling of subproject dependencies macro( "crest_find_package" @@ -24,12 +27,17 @@ macro( string(TOLOWER "${package}" _pkg_lc) string(TOUPPER "${package}" _pkg_uc) + + # iterate through lookup types in order foreach(method ${methods}) if(TARGET "${package}::${package}") break() endif() +######################################################################################### + + # look for a -config.cmake on the system if("${method}" STREQUAL "cmake") if(DEFINED "${_pkg_uc}_DIR") set("_${_pkg_uc}_DIR") @@ -42,6 +50,9 @@ macro( endif() endif() +######################################################################################### + + # look for dependency via pkgconf if("${method}" STREQUAL "pkgconf") find_package(PkgConfig QUIET) pkg_check_modules("${_pkg_uc}" QUIET "${package}") @@ -63,6 +74,9 @@ macro( endif() endif() +######################################################################################### + + # look for SOURCE in the subprojects directory (we usually prefer this one) if("${method}" STREQUAL "subproject") if(NOT DEFINED "${_pkg_uc}_SUBPROJECT") set("_${_pkg_uc}_SUBPROJECT") @@ -89,6 +103,9 @@ macro( endif() endif() +######################################################################################### + + # finally, we can try to download sources if("${method}" STREQUAL "fetch") message(STATUS "Retrieving ${package} from ${url}") include(FetchContent) @@ -112,6 +129,8 @@ macro( break() endif() +######################################################################################### + endforeach() if(TARGET "${package}::${package}") @@ -138,9 +157,10 @@ macro( endif() endmacro() -# +######################################################################################### +######################################################################################### + # Check current compiler version requirements. -# function (check_minimal_compiler_version lang compiler_versions) while(compiler_versions) list(POP_FRONT compiler_versions compiler version) diff --git a/docs/man/crest.adoc b/docs/man/crest.adoc index 7710aecb..b59a4d1b 100644 --- a/docs/man/crest.adoc +++ b/docs/man/crest.adoc @@ -275,6 +275,47 @@ The user is able to include additional constraints to **ALL** xtb**(1)** calcula Sum of population for structures considered in msRRHO average. [_default_: **0.9** (= 90%)] +=== options for MSREACT automated mass spectra fragment generator + + *-msreact*:: + start the msreact mode + + *-msnoattrh*:: + deactivate attractive potential between hydrogen and LMO centers + + *-msnshifts* _int_:: + perform n optimizations with randomly shifted atom postions (default 0) + + *-msnshifts* _int_:: + perform n optimizations with randomly shifted atom postions and repulsive potential applied to bonds (default 0) + + *-msnbonds* _int_:: + maximum number of bonds between atoms pairs for applying repulsive potential (default 3) + + *-msmolbar*:: + sort out topological duplicates by molbar codes (requires sourced "molbar") + + *-msinchi*:: + sort out topological duplicates by inchi codes (requires sourced "obabel") + + *-msnfrag* _int_:: + number of fragments that are printed by msreact (random selection) + + *-msiso*:: + print only non-dissociated structures (isomers) + + *-msnoiso*:: + print only dissociated structures + + *-mslargeprint*:: + do not remove temporary files and MSDIR with constrained optimizations + + *-chrg* _int_:: + set the molecules´ charge + + *-ewin* _float_:: + set energy window in for sorting out fragments kcal/mol, [default: 200.0 kcal/mol] + === Other tools for standalone use *-zsort*:: diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e407cfb5..51c70c1a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -27,6 +27,7 @@ add_subdirectory("rigidconf") add_subdirectory("discretize") add_subdirectory("entropy") add_subdirectory("legacy_algos") +add_subdirectory("msreact") set(dir "${CMAKE_CURRENT_SOURCE_DIR}") @@ -59,8 +60,6 @@ list(APPEND srcs "${dir}/marqfit.f90" "${dir}/minitools.f90" "${dir}/miscdata.f90" - "${dir}/msmod.f90" - "${dir}/msreact.f90" "${dir}/ncigeo.f90" "${dir}/ompmklset.F90" "${dir}/printouts.f90" diff --git a/src/classes.f90 b/src/classes.f90 index 800bc966..b6661674 100644 --- a/src/classes.f90 +++ b/src/classes.f90 @@ -471,6 +471,19 @@ module crest_data type(lwoniom_input),allocatable :: ONIOM_input !================================================! + !>--- msreact mode settings + logical :: msnoiso =.false. ! print only dissociated structures in msreact + logical :: msiso =.false. ! only print non-dissociated structures in msreact + logical :: msmolbar =.false. ! sort out duplicates by molbar + logical :: msinchi =.false. ! sort out duplicates by inchi + logical :: mslargeprint=.false. ! dont remove temporary files + logical :: msattrh=.true. ! add attractive potential for H-atoms + integer :: msnbonds = 3 ! distance of bonds up to nonds are stretched + integer :: msnshifts = 0 ! number of random shifts applied to whole mol + integer :: msnshifts2 = 0 ! number of random shifts applied to whole mol + integer :: msnfrag = 0 !number of fragments that are printed in msreact mode + character(len=80) :: msinput = '' !name of input file for msreact + !>--- general logical data logical :: allrot = .true. !> use all rotational constants for check instead of mean? logical :: altopt = .false. @@ -567,7 +580,6 @@ module crest_data logical :: water = .false. !> true if water is used as solvent (only QCG) logical :: wallsetup = .false. !> set up a wall potential? logical :: wbotopo = .false. !> set up topo with WBOs - contains procedure :: allocate => allocate_metadyn procedure :: deallocate => deallocate_metadyn diff --git a/src/confparse.f90 b/src/confparse.f90 index 7baa89b9..3538154c 100644 --- a/src/confparse.f90 +++ b/src/confparse.f90 @@ -305,6 +305,17 @@ subroutine parseflags(env,arg,nra) env%solv_file = '' env%solu_file = '' +!>--- options for msreact + env%msiso = .false. ! msiso and msnoiso are mutually exclusive !!! + env%msnoiso = .false. + env%msmolbar = .false. + env%mslargeprint=.false. ! dont remove temporary files + env%msattrh=.true. ! add attractive potential for H-atoms + env%msnbonds = 3 ! distance of bonds up to nonds are stretched + env%msnshifts = 0 ! number of random shifts applied to whole mol + env%msnshifts2 = 0 ! number of random shifts applied to whole mol + env%msnfrag = 0 !number of fragments that are printed in msreact mode + !&> !=========================================================================================! @@ -536,10 +547,11 @@ subroutine parseflags(env,arg,nra) env%autozsort = .false. exit - case ('-msreact') + case ('-msreact') env%crestver = crest_msreac env%preopt = .false. env%presp = .true. + env%ewin = 200.0d0 !> 200 kcal for msreact case ('-splitfile') ctmp = trim(arg(i+1)) @@ -1007,6 +1019,43 @@ subroutine parseflags(env,arg,nra) env%shake = nint(xx(1)) end select !> QCG end if + +!========================================================================================! +!------- Flags for msreact +!========================================================================================! + if (env%crestver == crest_msreac) then + select case (argument) !> msreact + case('-msnoiso') !> filter out non fragmentated structures in msreact + env%msnoiso=.true. + case('-msiso') !> filter out fragmentated structures in msreact + env%msiso=.true. + case('-msnbonds') ! give number of bonds up to which bias potential is added between atoms default 3 + call readl(arg(i + 1),xx,j) + env%msnbonds = xx(1) + case('-msnshifts') ! give number of times atoms are randomly shifted before optimization + call readl(arg(i + 1),xx,j) + env%msnshifts = xx(1) + case('-msnshifts2') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0 + call readl(arg(i + 1),xx,j) + env%msnshifts2 = xx(1) + case('-msnfrag') ! give number of structures that should be generated + call readl(arg(i + 1),xx,j) + env%msnfrag = xx(1) + case('-msmolbar') !> filter out structures with same molbar code in msreact + env%msmolbar=.true. + case('-msinchi') !> filter out structures with same inchi code in msreact + env%msinchi=.true. + case('-msnoattrh') !> add attractive potential for H-atoms + env%msattrh=.false. + case('-mslargeprint') !> additional printouts and keep MSDIR + env%mslargeprint=.true. + case('-msinput') ! give number of times atoms are randomly shifted before applying the constrained optimization default 0 + ctmp = trim(arg(i+1)) + if (ctmp(1:1) .ne. '-') then + env%msinput = trim(ctmp) + end if + end select !> msreact + end if !========================================================================================! !------- other general flags !========================================================================================! diff --git a/src/cregen.f90 b/src/cregen.f90 index 60683d4c..b53ec795 100644 --- a/src/cregen.f90 +++ b/src/cregen.f90 @@ -28,6 +28,7 @@ !> 6 - energy sorting only !> 9 - no sorting, only check groups !> 12 - no topology check, turn ewin to infty +!> 13 - no topology check, ewin and rmsd checking (msreact settings) !=========================================================================================! !=========================================================================================! @@ -364,6 +365,11 @@ subroutine cregen_files(env,fname,oname,cname,simpleset,userinput,iounit) oname = "crest_mecp_search.xyz.sorted" cname = "crest_ensemble.xyz" end if + if (simpleset == 13) then !> MSREACT files + fname = "crest_unique_products.xyz" + oname = "crest_unique_products.sorted" + cname = "crest_msreact_products.xyz" + end if if (simpleset == 15) then !> crossing files call checkname_xyz('confcross',fname,oname) cname = trim(fname)//'.unique' @@ -420,6 +426,13 @@ subroutine cregen_prout(env,simpleset,pr1,pr2,pr3,pr4) pr4 = .true. end if + if (simpleset == 13) then + pr1 = .false. + pr2 = .false. + pr3 = .false. + pr4 = .false. + end if + return end subroutine cregen_prout @@ -515,6 +528,21 @@ subroutine cregen_director(env,simpleset,checkbroken,sortE,sortRMSD,sortRMSD2, & anal = .false. end if + if(simpleset == 13)then !msreact mode + checkbroken = .false. + sorte = .true. + sortRMSD = .true. + sortRMSD2 = .false. + repairord = .false. + newfile = .true. + conffile = .true. + topocheck = .false. + checkez = .false. + bonusfiles= .false. + anal = .false. + + endif + return end subroutine cregen_director diff --git a/src/legacy_algos/cregen_old.f90 b/src/legacy_algos/cregen_old.f90 index 28208ff2..053f2ca4 100644 --- a/src/legacy_algos/cregen_old.f90 +++ b/src/legacy_algos/cregen_old.f90 @@ -32,7 +32,7 @@ subroutine cregen2(env) use axis_module use miscdata, only: rcov use utilities - !$ use omp_lib + use omp_lib implicit none type(systemdata) :: env ! MAIN STORAGE OS SYSTEM DATA diff --git a/src/legacy_algos/protonate.f90 b/src/legacy_algos/protonate.f90 index 1956ea9f..aec27abd 100644 --- a/src/legacy_algos/protonate.f90 +++ b/src/legacy_algos/protonate.f90 @@ -63,6 +63,14 @@ subroutine protonate(env,tim) logical :: ex + interface + subroutine xtblmo(env, print) + import :: systemdata + type(systemdata) :: env + logical, optional :: print + end subroutine xtblmo + end interface + !--- printout & clean directory call protclean call prothead @@ -87,7 +95,7 @@ subroutine protonate(env,tim) !--- do the xTB calculation for the LMOs call tim%start(1,'LMO calc.') - call xtblmo(env) + call xtblmo(env,.true.) call tim%stop(1) inquire (file='coordprot.0',exist=ex) if (.not.ex) then @@ -197,7 +205,7 @@ end subroutine protonate !-------------------------------------------------------------------------------------------- ! A quick single point xtb calculation and calculate LMOs !-------------------------------------------------------------------------------------------- -subroutine xtblmo(env) +subroutine xtblmo(env,print) use crest_parameters use iomod use crest_data @@ -207,6 +215,7 @@ subroutine xtblmo(env) character(len=:),allocatable :: jobcall integer :: io character(len=*),parameter :: pipe = ' > xtb.out 2>/dev/null' + logical, optional :: print ! leave the xtb.out file (e.g. for msreact mode) !---- setting threads if (env%autothreads) then @@ -228,7 +237,7 @@ subroutine xtblmo(env) !---- cleanup call remove(fname) - call remove('xtb.out') + if (.not. print) call remove('xtb.out') call remove('energy') call remove('charges') call remove('xtbrestart') diff --git a/src/meson.build b/src/meson.build index 09b3c439..e58da187 100644 --- a/src/meson.build +++ b/src/meson.build @@ -26,6 +26,7 @@ subdir('rigidconf') subdir('discretize') subdir('entropy') subdir('legacy_algos') +subdir('msreact') srcs += files( @@ -56,8 +57,6 @@ srcs += files( 'marqfit.f90', 'minitools.f90', 'miscdata.f90', - 'msmod.f90', - 'msreact.f90', 'ncigeo.f90', 'ompmklset.F90', 'printouts.f90', diff --git a/src/miscdata.f90 b/src/miscdata.f90 index 2b7ca31d..2c1a1866 100644 --- a/src/miscdata.f90 +++ b/src/miscdata.f90 @@ -69,6 +69,37 @@ module miscdata & * aatoau * 4.0_wp / 3.0_wp !&> +!&< + ! Radius used in QCxMS (in au) + real(wp), parameter :: QCxMS_Rad(118) = aatoau * [ & + & 0.32_wp,0.37_wp, & ! H,He + & 1.30_wp,0.99_wp,0.84_wp,0.75_wp,0.71_wp,0.64_wp,0.60_wp,0.62_wp, & ! Li-Ne + & 1.60_wp,1.40_wp,1.24_wp,1.14_wp,1.09_wp,1.04_wp,1.00_wp,1.01_wp, & ! Na-Ar + & 2.00_wp,1.74_wp, & ! K,Ca + & 1.59_wp,1.48_wp,1.44_wp,1.30_wp,1.29_wp, & ! Sc- + & 1.24_wp,1.18_wp,1.17_wp,1.22_wp,1.20_wp, & ! -Zn + & 1.23_wp,1.20_wp,1.20_wp,1.18_wp,1.17_wp,1.16_wp, & ! Ga-Kr + & 2.15_wp,1.90_wp, & ! Rb,Sr + & 1.76_wp,1.64_wp,1.56_wp,1.46_wp,1.38_wp, & ! Y- + & 1.36_wp,1.34_wp,1.30_wp,1.36_wp,1.40_wp, & ! -Cd + & 1.42_wp,1.40_wp,1.40_wp,1.37_wp,1.36_wp,1.36_wp, & ! In-Xe + & 2.38_wp,2.06_wp, & ! Cs,Ba + & 1.94_wp,1.84_wp,1.90_wp,1.88_wp,1.86_wp,1.85_wp,1.83_wp, & ! La-Eu + & 1.82_wp,1.81_wp,1.80_wp,1.79_wp,1.77_wp,1.77_wp,1.78_wp, & ! Gd-Yb + & 1.74_wp,1.64_wp,1.58_wp,1.50_wp,1.41_wp, & ! Lu- + & 1.36_wp,1.32_wp,1.30_wp,1.30_wp,1.32_wp, & ! -Hg + & 1.44_wp,1.45_wp,1.50_wp,1.42_wp,1.48_wp,1.46_wp, & ! Tl-Rn + & 2.42_wp,2.11_wp, & ! Fr,Ra + & 2.01_wp,1.90_wp,1.84_wp,1.83_wp,1.80_wp,1.80_wp,& ! Ac-Pu + ! from covalent 2009 covalent radii, such that it is complete up to 118 + & 1.49_wp, & ! Am + & 1.49_wp,1.51_wp,1.51_wp,1.48_wp,1.50_wp,1.56_wp,1.58_wp, & ! Cm-No + & 1.45_wp,1.41_wp,1.34_wp,1.29_wp,1.27_wp, & ! Lr- + & 1.21_wp,1.16_wp,1.15_wp,1.09_wp,1.22_wp, & ! -Cn + & 1.36_wp,1.43_wp,1.46_wp,1.58_wp,1.48_wp,1.57_wp ] ! Nh-Og +!&> + + !&< !> D3 pairwise van-der-Waals radii (only homoatomic pairs present here) real(wp), parameter :: vdw_D3(1:94) = aatoau * [& diff --git a/src/msmod.f90 b/src/msmod.f90 deleted file mode 100644 index 0a741de9..00000000 --- a/src/msmod.f90 +++ /dev/null @@ -1,173 +0,0 @@ -!================================================================================! -! This file is part of crest. -! -! Copyright (C) 2020 Philipp Pracht -! -! crest 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 3 of the License, or -! (at your option) any later version. -! -! crest 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 crest. If not, see . -!================================================================================! - -module msmod - - use iso_fortran_env, wp => real64 - implicit none - - - !-- storage for a single mol - type msmol - integer :: nat - integer,allocatable :: at(:) - real(wp),allocatable :: xyz(:,:) - real(wp) :: etot - integer :: chrg - integer :: nfrag - integer,allocatable :: frag(:) - integer :: gen - contains - procedure :: newmol => msmod_newmol - procedure :: dist => msmoldist - end type msmol - - public :: msmol - - !-- isomer list - type ilist - integer :: nmol - type(msmol),allocatable :: mol(:) - logical,allocatable :: new(:) - contains - procedure :: dealloc => deallocate_ilist - procedure :: append => ilist_append - end type ilist - - !-- fragmentized structure list - type flist - integer :: nmol - type(msmol),allocatable :: mol(:) - logical,allocatable :: new(:) - contains - procedure :: dealloc => deallocate_flist - end type flist - - - !-- global data for msreact - type msobj - - real(wp) :: ethr = 1.0_wp !cov. struc. chekc. thr. - real(wp) :: wbormsd = 0.5_wp !wbo rmsd comp. thr. - - real(wp) :: T = 3000.0_wp !start temperature - real(wp) :: fc = 0.05_wp !start fc - real(wp) :: cdist = 1.5_wp !constraing distance scaling factor of rcov - integer :: maxc = 30 !max optimization cycle - - type(ilist) :: il - type(flist) :: fl - - end type msobj - - public :: msobj - - private -contains - -subroutine deallocate_ilist(self) - class(ilist) :: self - if(allocated(self%mol)) deallocate(self%mol) - if(allocated(self%new)) deallocate(self%new) - return -end subroutine deallocate_ilist - -subroutine deallocate_flist(self) - class(flist) :: self - if(allocated(self%mol)) deallocate(self%mol) - if(allocated(self%new)) deallocate(self%new) - return -end subroutine deallocate_flist - -subroutine msmod_newmol(self,nat,xyz,at,etot,chrg,gen) - implicit none - class(msmol) :: self - integer :: nat - real(wp) :: xyz(3,nat) - integer :: at(nat) - real(wp) :: etot - integer :: chrg - integer :: gen - self%nat = nat - allocate(self%xyz(3,nat)) - self%xyz = xyz - allocate(self%at(nat)) - self%at = at - self%etot = etot - self%chrg = chrg - self%gen = gen - return -end subroutine msmod_newmol - -!calculate the distance between two atoms i and j -function msmoldist(self,i,j) result(dist) - implicit none - class(msmol) :: self - integer :: i,j - real(wp) :: dist - dist = 0.0_wp - if(allocated(self%xyz))then - dist=(self%xyz(1,i)-self%xyz(1,j))**2 + & - & (self%xyz(2,i)-self%xyz(2,j))**2 + & - & (self%xyz(3,i)-self%xyz(3,j))**2 - dist = sqrt(dist) - endif - return -end function msmoldist - - - -subroutine ilist_append(self,nat,at,xyz,etot,chrg,gen) - implicit none - class(ilist) :: self - type(msmol) :: mol - integer :: nat - real(wp) :: xyz(3,nat) - integer :: at(nat) - real(wp) :: etot - integer :: chrg - integer :: gen - - type(msmol),allocatable :: dummy(:) - logical,allocatable :: btmp(:) - integer :: n - call mol%newmol(nat,xyz,at,etot,chrg,gen) - if(.not.allocated(self%mol))then - self%nmol = 1 - allocate(self%mol(1)) - allocate(self%new(1)) - self%mol(1) = mol - self%new(1) = .true. - else - n = self%nmol + 1 - allocate(dummy(n)) - allocate(btmp(n)) - dummy(1:self%nmol) = self%mol(1:self%nmol) - dummy(n) = mol - btmp(1:self%nmol) = self%new(1:self%nmol) - btmp(n) = .true. - self%nmol = n - call move_alloc(btmp, self%new) - call move_alloc(dummy, self%mol) - endif - return -end subroutine ilist_append - - -end module msmod diff --git a/src/msreact.f90 b/src/msreact.f90 deleted file mode 100644 index dc775369..00000000 --- a/src/msreact.f90 +++ /dev/null @@ -1,346 +0,0 @@ -!================================================================================! -! This file is part of crest. -! -! Copyright (C) 2020 Stefan Grimme, Philipp Pracht -! -! crest 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 3 of the License, or -! (at your option) any later version. -! -! crest 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 crest. If not, see . -!================================================================================! - -!==============================================================! -! the handler function is the sub-program that is called -! by crest, i.e., within this routine the seperate input file is -! read and the calculation is started. -!==============================================================! -subroutine msreact_handler(env,tim) - use crest_parameters - use crest_data - use msmod - use strucrd - implicit none - - type(systemdata) :: env - type(timer) :: tim - - type(msobj) :: mso !a new msreact object - type(coord) :: struc - - interface - subroutine msreact_topowrap(mol,pair,paths,wboname) - import :: msmol - implicit none - type(msmol) :: mol - integer :: pair(mol%nat*(mol%nat+1)/2) - integer :: paths(mol%nat*(mol%nat+1)/2,mol%nat) - character(len=*),optional :: wboname - end subroutine msreact_topowrap - end interface - integer :: nat - integer,allocatable :: pair(:) - integer,allocatable :: paths(:,:) - integer :: k - - call tim%start(1,'MSREACT') - - !-- read the input coord and put it into the - ! iso-list as Gen 0 structure - call struc%open('coord') - struc%xyz=struc%xyz*bohr !to Angström, from this point on by convention! - call mso%il%append(struc%nat,struc%at,struc%xyz,0.0_wp,env%chrg,0) - call struc%deallocate() - - !-- additional input file could be read here - !call msinputreader(mso) - - - nat = mso%il%mol(1)%nat - k=nat*(nat+1)/2 - allocate(pair(k),paths(k,nat)) - call msreact_topowrap(mso%il%mol(1),pair,paths,'wbo') - - !-- setting the threads for correct parallelization - if(env%autothreads)then - call ompautoset(env%threads,6,env%omp,env%MAXRUN,0) !set global OMP/MKL variable for xtb jobs - endif - - !-- do the directory setup and optimizations - call msreact(mso,mso%il%mol(1),nat,pair,3) - - - deallocate(paths,pair) - - call tim%stop(1) - return -end subroutine msreact_handler - -!==============================================================! -! the main implementation of the msreact algo should go here -!==============================================================! -subroutine msreact(mso,mol,nat,pair,nbonds) - use crest_parameters - use msmod - use iomod - use miscdata, only: rcov - use utilities, only: lin - implicit none - - type(msobj) :: mso !main storage object - type(msmol) :: mol ! xyz etc - - integer :: nat - integer :: pair(nat*(nat+1)/2) - !integer :: paths(nat*(nat+1)/2,nat) - integer :: nbonds - integer :: i,j,k - integer :: p - integer :: np - integer :: io - - character(len=:),allocatable :: subdir - character(len=40) :: pdir - character(len=512) :: thisdir - real(wp) :: constr_dist - - !-- main subdirectory handling - call getcwd(thisdir) - subdir='MSDIR' - io = makedir(subdir) - call chdir(subdir) - - - !-- get specific pairs - np=0 - do p=1,nbonds ! bonds in between -! write(*,'(1x,a,i0,a)') '1,',p+1,' pairs' - do i=1,nat - do j=i,nat - k=lin(i,j) - if(p.eq.1.and.pair(k).eq.1.and.(mol%at(i).eq.1.or.mol%at(j).eq.1)) cycle - if(pair(k)==p)then - np = np+1 - write(pdir,'(a,i0)')'Pair_',np - constr_dist = mso%cdist*(rcov(mol%at(i))+rcov(mol%at(j)))*bohr + float(p) -! write(*,*) mol%at(i),mol%at(j),constr_dist - call isodir(mso,trim(pdir),mol,i,j,constr_dist) - endif - enddo - enddo - enddo - - write(*,*) '# of distortions',np - call msreact_jobber(np,'Pair_',.false.) - - call msreact_collect(mol%nat,np,'products.xyz') - call rename(subdir//'/'//'products.xyz','products.xyz') - call chdir(thisdir) - return -end subroutine msreact - - -!============================================================! -! make a dir for a structure without fragments, -! a controlfile with constraints on atoms A and B (at dist D) -! will be written into the directory -!============================================================! -subroutine isodir(mso,dirname,mol,A,B,D) - use crest_parameters - use msmod - use iomod - use strucrd, only : wrxyz - implicit none - type(msobj) :: mso - character(len=*) :: dirname - type(msmol) :: mol - integer :: A,B - real(wp) :: D - - character(len=:),allocatable :: fname - character(len=20) :: dumm - integer :: io,ich - - io = makedir(dirname) !create the directory - - fname = trim(dirname)//'/'//'struc.xyz' - open(newunit=ich,file=fname) - call wrxyz(ich,mol%nat,mol%at,mol%xyz) - close(ich) - - fname = trim(dirname)//'/'//'.CHRG' - open(newunit=ich,file=fname) - write(ich,'(i0)') mol%chrg + 1 ! EI +1, DEA -1, CID 0 - close(ich) - - fname = trim(dirname)//'/'//'.xc1' - open(newunit=ich, file=fname) - write(ich,'(a)') '$scc' - write(dumm,'(f16.2)') mso%T - write(ich,'(1x,a,a)')'temp=',adjustl(trim(dumm)) - write(ich,'(a)')'$constrain' - write(dumm,'(f16.4)') mso%fc - write(ich,'(3x,a,a)')'force constant=',adjustl(trim(dumm)) - write(ich,'(3x,a,1x,i0,a,1x,i0,a,1x,f8.5)') 'distance:',A,',',B,',',D - close(ich) - - fname = trim(dirname)//'/'//'.xc2' - open(newunit=ich, file=fname) - write(ich,'(a)') '$scc' - write(dumm,'(f16.2)') mso%T - write(ich,'(1x,a,a)')'temp=',adjustl(trim(dumm)) - write(ich,'(a)') '$opt' - write(ich,'(1x,a)') 'maxcycle=5' - write(ich,'(a)') '$write' - write(ich,'(1x,a)') 'wiberg=true' - close(ich) - - return -end subroutine isodir - -!=====================================================================! -! The job construction routine for MSREACT -! (will have to be modified later, for now it is for testing) -!=====================================================================! -subroutine msreact_jobber(ndirs,base,niceprint) - use crest_parameters - use msmod - use iomod - implicit none - integer :: ndirs - character(len=*) :: base - logical :: niceprint - - character(len=1024) :: jobcall - character(len=1024) :: jobcall2 - - jobcall = '' - jobcall2 = '' - - write(jobcall,'(a)') 'xtb struc.xyz --opt loose --input .xc1 > split.out 2>/dev/null' - write(jobcall2,'(a)') 'xtb xtbopt.xyz --opt crude --input .xc2 > xtb.out 2>/dev/null' - jobcall = trim(jobcall)//' ; '//trim(jobcall2) - - !-- directories must be numbered consecutively - call opt_OMP_loop(ndirs,base,jobcall,niceprint) - write(*,*) - write(*,*) 'done.' - return -end subroutine msreact_jobber - - - -!=====================================================================! -! A wrapper to generate the topology for a molecule within the -! MSREACT subprogram -!=====================================================================! -subroutine msreact_topowrap(mol,pair,paths,wboname) - use crest_parameters - use msmod - use zdata - use adjacency - use utilities, only: lin - implicit none - type(msmol) :: mol - integer :: pair(mol%nat*(mol%nat+1)/2) - !integer :: pair(mol%nat,mol%nat) - integer :: paths(mol%nat*(mol%nat+1)/2,mol%nat) - character(len=*),optional :: wboname - type(zmolecule) :: zmol - - integer,allocatable :: A(:,:) - integer,allocatable :: prev(:,:) - real(wp),allocatable :: E(:,:) - real(wp),allocatable :: dist(:,:) - - integer :: lpath,i,j,k - integer,allocatable :: path(:) - logical :: ex - - - ex=.false. - if(present(wboname))then - inquire(file=wboname,exist=ex) - endif - if(ex)then - call simpletopo(mol%nat,mol%at,mol%xyz,zmol,.false.,.false.,wboname) - else - mol%xyz = mol%xyz / bohr !CN based topo requires Bohrs - call simpletopo(mol%nat,mol%at,mol%xyz,zmol,.false.,.false.,'') - mol%xyz = mol%xyz * bohr - endif - - allocate(A(mol%nat,mol%nat),E(mol%nat,mol%nat)) - call zmol%adjacency(A,E) - - allocate(prev(mol%nat,mol%nat),dist(mol%nat,mol%nat)) - - call FloydWarshall(mol%nat,A,E,dist,prev) - allocate(path(mol%nat), source = 0) - do i=1,mol%nat - do j=i,mol%nat - path = 0 - call getPathFW(mol%nat,prev,i,j,path,lpath) - !write(*,*) path(1:lpath) - k=lin(i,j) - pair(k) = lpath - 1 ! number of bonds - paths(k,:) = path(:) - enddo - enddo - - deallocate(dist,prev) - deallocate(E,A) - - call zmol%deallocate() !clear the zmol memory - return -end subroutine msreact_topowrap - -!========================================================================! -! collect structures of optimized molecules -! xyz files should still have the same number and order of atoms -!========================================================================! -subroutine msreact_collect(nat,np,outfile) - use crest_parameters - use strucrd - implicit none - integer :: nat - integer :: np - character(len=*) :: outfile - integer :: ich - character(len=40) :: pdir - character(len=:),allocatable :: optfile - character(len=128) :: newcomment - integer :: p,p2 - logical :: ex - integer,allocatable :: at(:) - real(wp),allocatable :: xyz(:,:) - real(wp) :: etot - - - allocate(at(nat),xyz(3,nat)) - open(newunit=ich,file=outfile) - p=0 - do p2=1,np - write(pdir,'(i0,i0,a,i0)')1,p+1,'Pair_',p2 - write(pdir,'(a,i0)')'Pair_',p2 - optfile=trim(pdir)//'/'//'xtbopt.xyz' - inquire(file=optfile,exist=ex) - if(ex)then - call rdcoord(optfile,nat,at,xyz,etot) - xyz = xyz*bohr - write(newcomment,'(1x,f18.8,5x,a)')etot,trim(pdir) - call wrxyz(ich,nat,at,xyz,newcomment) - endif - enddo - close(ich) - - deallocate(xyz,at) - return -end subroutine msreact_collect diff --git a/src/msreact/CMakeLists.txt b/src/msreact/CMakeLists.txt new file mode 100644 index 00000000..b65ec325 --- /dev/null +++ b/src/msreact/CMakeLists.txt @@ -0,0 +1,24 @@ +# This file is part of crest. +# SPDX-Identifier: LGPL-3.0-or-later +# +# crest 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 3 of the License, or +# (at your option) any later version. +# +# crest 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 crest. If not, see . + +set(dir "${CMAKE_CURRENT_SOURCE_DIR}") + +list(APPEND srcs + "${dir}/msreact.f90" + "${dir}/msmod.f90" +) + +set(srcs ${srcs} PARENT_SCOPE) diff --git a/src/msreact/meson.build b/src/msreact/meson.build new file mode 100644 index 00000000..4f3cafe6 --- /dev/null +++ b/src/msreact/meson.build @@ -0,0 +1,20 @@ +# This file is part of crest. +# SPDX-Identifier: LGPL-3.0-or-later +# +# crest 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 3 of the License, or +# (at your option) any later version. +# +# crest 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 crest. If not, see . + +srcs += files( + 'msreact.f90', + 'msmod.f90', +) diff --git a/src/msreact/msmod.f90 b/src/msreact/msmod.f90 new file mode 100644 index 00000000..ba25c1e0 --- /dev/null +++ b/src/msreact/msmod.f90 @@ -0,0 +1,1154 @@ +!================================================================================! +! This file is part of crest. +! +! Copyright (C) 2020-2024 Philipp Pracht, Johannes Gorges +! +! crest 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 3 of the License, or +! (at your option) any later version. +! +! crest 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 crest. If not, see . +!================================================================================! + +module msmod + use crest_parameters + use crest_data + use strucrd + use iomod + implicit none + + !>-- storage for a single mol + type msmol + integer :: nat + integer,allocatable :: at(:) + real(wp),allocatable :: xyz(:,:) + real(wp) :: etot + integer :: chrg + integer :: nfrag + !> index of fragatoms 1 first fragment 2 second fragment + integer,allocatable :: fragi(:) + !> number of fragments in product + !> 1: isomer 1< : fragmentation; 0, not determined yet, -1 means we sort this out + integer :: fragcount + real(wp) :: molmass + contains + procedure :: newmol => msmod_newmol + procedure :: dist => msmoldist + end type msmol + + public :: msmol + + !-- product list + type plist + integer :: nmol + type(msmol),allocatable :: mol(:) + logical,allocatable :: new(:) + contains + procedure :: dealloc => deallocate_plist + procedure :: append => plist_append + procedure :: remove => plist_remove + end type plist + + !>-- global data for msreact + type msobj + real(wp) :: rcut = 1.3_wp !> cutoff for factor of sum of covalent radii for detection of fragmentation of fragment_structures + real(wp) :: T = 5000.0_wp !> electronic temperature! favours open-shell fragment generation + real(wp) :: fc = 0.05_wp !> repulsive force constant in constrained optimizations + real(wp) :: fc_attr = -0.05_wp !> attractive force constant in constrained optimizations + real(wp) :: cdist = 1.5_wp !> constraing distance scaling factor of rcov + real(wp) :: cdist_att = 0.5_wp !> constraing distance scaling factor of rcov for attractive part + real(wp) :: distthr_att = 4.0_wp !> distance threshold in Angstrom for attractive part of constraint + real(wp) :: atomshift = 0.75_wp !> shift of atoms in random atom displacement + real(wp) :: fragdist = 0.0_wp !> increase distance by to separate fragments (useful for subsequent transition state search) + integer :: maxc = 15 !>max optimization cycle in constrained optimizations + + type(plist) :: pl + end type msobj + + public :: msobj + +!=======================================================================================! +!=======================================================================================! +contains !> MODULE PROCEDURES START HERE +!=======================================================================================! +!=======================================================================================! + + subroutine deallocate_plist(self) + class(plist) :: self + if (allocated(self%mol)) deallocate (self%mol) + if (allocated(self%new)) deallocate (self%new) + return + end subroutine deallocate_plist + +!=======================================================================================! + + subroutine msmod_newmol(self,nat,xyz,at,etot,chrg,fragcount,molmass) + implicit none + class(msmol) :: self + integer :: nat + real(wp) :: xyz(3,nat) + integer :: at(nat) + real(wp) :: etot + integer :: chrg + integer :: fragcount + real(wp) :: molmass + self%nat = nat + allocate (self%xyz(3,nat)) + self%xyz = xyz + allocate (self%at(nat)) + self%at = at + self%etot = etot + self%chrg = chrg + self%fragcount = fragcount + self%molmass = molmass + return + end subroutine msmod_newmol + +!=======================================================================================! + + function msmoldist(self,i,j) result(dist) +!*************************************************** +!* calculate the distance between two atoms i and j +!*************************************************** + implicit none + class(msmol) :: self + integer :: i,j + real(wp) :: dist + dist = 0.0_wp + if (allocated(self%xyz)) then + dist = (self%xyz(1,i)-self%xyz(1,j))**2+ & + & (self%xyz(2,i)-self%xyz(2,j))**2+ & + & (self%xyz(3,i)-self%xyz(3,j))**2 + dist = sqrt(dist) + end if + return + end function msmoldist + +!=======================================================================================! + + subroutine plist_append(self,nat,at,xyz,etot,chrg,fragcount,molmass) + implicit none + class(plist) :: self + type(msmol) :: mol + integer :: nat + real(wp) :: xyz(3,nat) + integer :: at(nat) + real(wp) :: etot + integer :: chrg + integer :: fragcount + real(wp) :: molmass + + type(msmol),allocatable :: dummy(:) + logical,allocatable :: btmp(:) + integer :: n + call mol%newmol(nat,xyz,at,etot,chrg,fragcount,molmass) + if (.not.allocated(self%mol)) then + self%nmol = 1 + allocate (self%mol(1)) + allocate (self%new(1)) + self%mol(1) = mol + self%new(1) = .true. + else + n = self%nmol+1 + allocate (dummy(n)) + allocate (btmp(n)) + dummy(1:self%nmol) = self%mol(1:self%nmol) + dummy(n) = mol + btmp(1:self%nmol) = self%new(1:self%nmol) + btmp(n) = .true. + self%nmol = n + call move_alloc(btmp,self%new) + call move_alloc(dummy,self%mol) + end if + return + end subroutine plist_append + +!=======================================================================================! + + subroutine plist_remove(self,index) +!**************************************************** +!* remove entry from product with index i from list +!**************************************************** + implicit none + class(plist) :: self + type(msmol) :: mol + + integer :: n,i + integer,intent(in) :: index + + if (index < 1.or.index > self%nmol) then + write (*,*) "Invalid index. Cannot remove molecule." + return + end if + + do i = index,self%nmol-1 + self%mol(i) = self%mol(i+1) ! replace the entry with the next one + end do + + self%nmol = self%nmol-1 ! ignore last entry + + end subroutine plist_remove + +!=======================================================================================! + + subroutine get_input_energy(env,etemp,etot) +!****************************************************************************** +!* A quick xtb geometry optimization in xyz coordinates to get starting energy +!****************************************************************************** + implicit none + type(systemdata) :: env + character(len=80) :: fname,pipe + character(len=:),allocatable :: jobcall + logical :: fin + character(len=256) :: atmp + integer :: ich,iost,io,i + type(coord) :: mol + integer :: ntopo + integer,allocatable :: topo(:) + real(wp),intent(out) :: etot + real(wp) :: etemp ! electronic temperature in K + logical :: tchange = .false. + logical :: ldum + +!---- setting threads + if (env%autothreads) then + call ompautoset(env%threads,7,env%omp,env%MAXRUN,1) !set the global OMP/MKL variables for the xtb jobs + end if + +!---- small header + write (*,*) + call smallhead('xTB Geometry Optimization') +!---- some options + pipe = ' > xtb.out 2>/dev/null' + call remove('gfnff_topo') + if (.not.env%chargesfile) call remove('charges') + call remove('grad') + call remove('mos') + call remove('xtbopt.log') + call remove('xtbrestart') + + fname = '.CHRG' + open (newunit=ich,file=fname) + write (ich,'(i0)') env%chrg ! EI +1, DEA -1, CID 0, give in crest call + close (ich) + +!---- input xyz file + fname = env%inputcoords + +! write (jobcall,'(a,1x,a,f10.4,1x,a,1x,a)') & +! & trim(env%ProgName),trim(fname)//" --sp --etemp ",etemp,trim(env%gfnver),trim(pipe) + jobcall = trim(env%ProgName) + jobcall = trim(jobcall)//' '//trim(fname) + write(atmp,'(f10.4)') etemp + jobcall = trim(jobcall)//' --sp --etemp '//trim(atmp) + jobcall = trim(jobcall)//trim(env%gfnver)//trim(pipe) + + call execute_command_line(trim(jobcall),exitstat=io) + + call minigrep('xtb.out','finished run',fin) + if (.not.fin) then + write (*,*) + write (*,*) ' Initial singlepoint calculation failed!' + write (*,*) ' Please check your input.' + error stop + end if + write (*,*) 'Input energy successfully calculated.' + + call grepval('xtb.out',"| TOTAL ENERGY",ldum,etot) +!---- cleanup + call remove('xtb.out') + call remove('energy') + if (.not.env%chargesfile) call remove('charges') + call remove('grad') + call remove('mos') + call remove('xtbopt.log') + call remove('xtbtopo.mol') + call remove('.xtbopttok') + call remove('xtbrestart') + call remove('gfnff_topo') + end subroutine get_input_energy + +!=======================================================================================! + + subroutine readbasicpos(env,nbaseat,basicatlist) +!******************************************************************** +!* this subroutine reads the lmo.out file of an xtb -lmo calculation +!* and identifies the atoms which have a LP or pi-Orbital +!******************************************************************** + implicit none + type(systemdata) :: env + integer :: i,io,ich,at1,at2,j + integer,intent(out) :: nbaseat + integer,intent(out),allocatable :: basicatlist(:) + integer,allocatable :: dumlist(:) + character(len=512) :: tmp + character(len=64) :: fname + character(len=64) :: type,dumc + real(wp) :: dumr + logical :: ex + fname = 'xtb.out' + inquire (file=fname,exist=ex) + if (.not.ex) then + write (*,*) "xtb.out not found" + stop + end if + nbaseat = 0 + open (newunit=ich,file=fname) + do + read (ich,'(a)',iostat=io) tmp + if (index(tmp,'files:') .ne. 0) exit + if (index(tmp,'pi ') .ne. 0) nbaseat = nbaseat+2 ! count first two atoms of pi or delpi bond with highest participation and ignore the rest + if (index(tmp,'LP ') .ne. 0) nbaseat = nbaseat+1 ! count LP + end do + allocate (dumlist(nbaseat)) + dumlist = 0 + j = 0 + rewind (ich) + do + read (ich,'(a)',iostat=io) tmp + if (index(tmp,'starting deloc pi') .ne. 0) exit + if (index(tmp,'files:') .ne. 0) exit + if (index(tmp,'pi ') .ne. 0) then + backspace (ich) + if (tmp(64:64) == ' ') then ! if element symbol has one character, there is a space and we use this routine + if (tmp(80:80) == ' ') then + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumc,dumr,at2,dumc,dumc,dumr ! first element has one and second has one character + else + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumc,dumr,at2,dumc,dumr ! first element has one and second has two character + end if + else + if (tmp(80:80) == ' ') then + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumr,at2,dumc,dumc,dumr ! first element has two and second has pm character + else + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumr,at2,dumc,dumr ! first element has two and second has two character + end if + end if + if (findloc(dumlist,at1,1) .eq. 0) then + j = j+1 + dumlist(j) = at1 + end if + if (findloc(dumlist,at2,1) .eq. 0) then + j = j+1 + dumlist(j) = at2 + end if + end if + + if (index(tmp,'LP ') .ne. 0) then + backspace (ich) + if (tmp(64:64) == ' ') then ! if element symbol has one character, there is a space and we use this routine + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumc,dumr + else ! if element symbol has two characters, there is no space and we use this routine + read (ich,*) dumr,type,dumr,dumr,dumr,dumr,dumr,at1,dumc,dumr + end if + if (findloc(dumlist,at1,1) .eq. 0) then + j = j+1 + dumlist(j) = at1 + end if + end if + end do + close (ich) + basicatlist = pack(dumlist,dumlist .ne. 0) ! sort out zeroes + nbaseat = size(basicatlist) + write (*,"('Protonation sites at atom positions:', *(i5))") (basicatlist(i),i=1,nbaseat) + call remove('xtb.out') + end subroutine readbasicpos + +!=======================================================================================! + + subroutine fragment_structure(nat,oz,xyz,rcut,at1,at2,frag,fragcount) +!***************************************************************************************** +!* Taken and slightly modified from QCxMS mass spectra code https://github.com/qcxms/QCxMS +!* Subroutine for definition of two or more fragments +!* if at1 = 0 : global separation in (nonbonded parts), beginning with atom at2 +!* if at1 and at2 <> 0 : define fragments only if a at1-at2 bond (and no other) exists +!* if at1 and at2 = 0 : delete all fragment assignments +!* no bond if rcut times the cov.dist. +!* works better than mrec +!***************************************************************************************** + use miscdata,only:Rad => QCxMS_Rad + implicit none + integer :: at1,at2,nat + integer :: i,j + integer :: attotal,currentfrag + integer :: oz(nat),frag(nat) + integer :: fragcount + real(wp),intent(in) :: xyz(3,nat) + real(wp) :: rcov,r + real(wp) :: rcut + + logical :: finish + logical,allocatable :: connect(:,:) + + allocate (connect(nat,nat)) + connect(1:nat,1:nat) = .false. + + do i = 1,nat-1 + do j = i+1,nat + r = sqrt((xyz(1,i)-xyz(1,j))**2+(xyz(2,i)-xyz(2,j))**2 & + & +(xyz(3,i)-xyz(3,j))**2) + rcov = rcut*0.5_wp*(Rad(oz(i))+Rad(oz(j))) + if (r .lt. rcov) then + connect(i,j) = .true. + connect(j,i) = .true. + end if + end do + end do + if ((at1 .eq. 0).and.(at2 .eq. 0)) then + do i = 1,nat + frag(i) = 1 + end do + return + else + + do i = 1,nat + frag(i) = 0 + end do + + frag(at1) = 1 + attotal = 1 + + if (at2 .ne. 0) then + connect(at1,at2) = .false. + connect(at2,at1) = .false. + end if + + finish = .false. + currentfrag = 0 + + do while (attotal .ne. nat) + + currentfrag = currentfrag+1 + + ! cycle through atoms and find connected ones + + do while (.not. (finish)) + finish = .true. + do i = 1,nat + if (frag(i) .eq. currentfrag) then + do j = 1,nat + if (connect(i,j)) then + if (frag(j) .eq. 0) then + frag(j) = currentfrag + attotal = attotal+1 + finish = .false. + elseif (frag(j) .eq. currentfrag) then + cycle + end if + end if + end do + end if + end do + end do + + ! find the first atom in the next fragment + + do i = 1,nat + if (frag(i) .eq. 0) then + frag(i) = currentfrag+1 + attotal = attotal+1 + exit + end if + end do + finish = .false. + end do + + end if + + do i = 1,3 ! is enough, we only need to know we have 1, 2 or more than 2 fragments + if (count(frag == i) .gt. 0) then + fragcount = i + end if + end do + + deallocate (connect) + return + end subroutine fragment_structure + +!=======================================================================================! + + subroutine sortoutduplicates(env,mso) +!****************************************************************************** +!* sort out duplicates with molbar or inchi by openbabel (needs to be sourced) +!****************************************************************************** + implicit none + type(msobj) :: mso + integer :: np + integer :: i,j + integer :: rm !removecounter + character(len=1024),allocatable :: barcodes(:) + character(len=2048) :: duplicates + logical,allocatable :: double(:) + type(systemdata) :: env ! MAIN STORAGE OS SYSTEM DATA + logical :: lprint + logical :: stat + + lprint = env%mslargeprint + np = mso%pl%nmol + allocate (double(np),barcodes(np)) + + ! compute barcodes + if (env%msmolbar) then + call calcmolbar(env,mso,np,barcodes,stat) + if (.not.stat) return + elseif (env%msinchi) then + call calcinchi(env,mso,np,barcodes,stat) + if (.not.stat) return + else + return ! no topology check for duplicated structures selected + end if + + ! identify duplicates + double = .false. + do i = 1,np + if (double(i)) cycle + duplicates = "" + do j = i+1,np + if (.not.double(j)) then + if (trim(barcodes(i)) == trim(barcodes(j))) then + double(j) = .true. + if (lprint) write (duplicates,'(a,i0)') trim(duplicates)//" ",j + else + cycle + end if + end if + end do + if (lprint.and.len(trim(duplicates)) .gt. 0) then + write (*,'(a,i0,a)') "removed duplicates of ",i,": "//trim(duplicates) + end if + end do + + ! remove duplicates + rm = 0 + do i = 1,np + rm = rm+1 + if (double(i)) then + call mso%pl%remove(rm) + rm = rm-1 + end if + end do + + write (*,*) "sorted out ",np-mso%pl%nmol," fragments accoding to topology check" + write (*,*) "number of products is now:",mso%pl%nmol-1 + + end subroutine sortoutduplicates + +!=======================================================================================! + + subroutine calcmolbar(env,mso,np,barcodes,stat) +!****************************************************** +!* compute molbar codes for all structures with molbar +!****************************************************** + implicit none + type(systemdata) :: env ! MAIN STORAGE OS SYSTEM DATA + type(msobj) :: mso + integer :: np + integer :: i,ich,iocheck,io + character(len=:),allocatable :: prog + character(len=:),allocatable :: subdir + character(len=1024) :: barcodes(np) + character(len=512) :: jobcall + character(len=512) :: thisdir + character(len=80) :: fname + logical :: sourced,stat + + stat = .true. + + prog = 'molbar' + call checkprog_silent(prog,.true.,io) + if (io .ne. 0) then + write (*,*) trim(prog)," not found, no topology check for duplicated structures possible" + stat = .false. + return + end if + + call getcwd(thisdir) + subdir = 'topodir' + io = makedir(subdir) + call chdir(subdir) + + do i = 1,np + write (fname,'(i0,a)') i,"temp.xyz" + open (newunit=ich,file=fname) + call wrxyz(ich,mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),mso%pl%mol(i)%etot) + close (ich) + end do + + ! molbar on multiple core is faster but was unstable in some cases, single core more stable + !write(jobcall,'(a)') 'molbar *temp.xyz -s > molbar.out' + + write (jobcall,'(a,i0,a)') 'molbar *temp.xyz -s -T ',env%threads,' > molbar.out' + write (*,*) "Calling molbar for sorting out duplicates by molbar" + + call execute_command_line(trim(jobcall)) + + do i = 1,np + write (fname,'(i0,a)') i,'temp.mb' + open (newunit=ich,file=fname,status="old",action="read") + read (ich,'(a)',iostat=iocheck) barcodes(i) + close (ich) + end do + + call chdir(thisdir) + call rmrf('topodir') + end subroutine calcmolbar + +!=======================================================================================! + + subroutine calcinchi(env,mso,np,barcodes,stat) +!******************************************************** +!* compute inchi codes for all structures with openbabel +!******************************************************** + implicit none + type(systemdata) :: env ! MAIN STORAGE OS SYSTEM DATA + type(msobj) :: mso + integer :: np + integer :: i,ich,iocheck,io + character(len=:),allocatable :: prog + character(len=:),allocatable :: subdir + character(len=1024) :: barcodes(np) + character(len=512) :: jobcall + character(len=512) :: thisdir + character(len=80) :: fname + logical :: sourced,stat + + stat = .true. + + prog = 'obabel' + call checkprog_silent(prog,.true.,io) + if (io .ne. 0) then + write (*,*) trim(prog)," not found, no topology check for duplicated structures possible" + stat = .false. + return + end if + + call getcwd(thisdir) + subdir = 'topodir' + io = makedir(subdir) + call chdir(subdir) + + open (newunit=ich,file='temp.xyz') + do i = 1,np + call wrxyz(ich,mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),mso%pl%mol(i)%etot) + end do + close (ich) + + write (jobcall,'(a)') 'obabel -i xyz temp.xyz -o inchi > inchicodes 2>obabelout' + write (*,*) "Calling obabel for sorting out duplicates by InChi codes" + call execute_command_line(trim(jobcall)) + + open (newunit=ich,file='inchicodes',status="old",action="read") + do i = 1,np + read (ich,"(a)") barcodes(i) + end do + close (ich) + + call chdir(thisdir) + call rmrf('topodir') + end subroutine calcinchi + +!=======================================================================================! + + subroutine write_fragments(env,mso,estart,nisomers,nfragpairs,fname,lprint) +!********************************************************* +!* write fragments to files and optionally to directories +!********************************************************* + use atmasses + implicit none + character(len=*) :: fname + real(wp),intent(in) :: estart + integer,intent(in) :: nisomers,nfragpairs + type(systemdata) :: env + type(msobj) :: mso + integer :: prstruc ! number of structures to print + integer :: incr + integer :: np0,np,np2 + integer :: i,j,k,ich1,ich2,ich,r + integer :: nat + integer :: natf(2) ! number of atoms of fragments 1 and 2 + integer,allocatable :: atf(:,:) ! atomtypes of fragment + character(len=1024) :: thisdir,isodir,fragdir + character(len=40) :: strucname + character(len=80) :: comment + logical :: ex + logical :: lprint + character(len=40) :: sumform,sumformula + real(wp) :: mass,erel + + nat = mso%pl%mol(1)%nat + call getcwd(thisdir) + np = mso%pl%nmol + np0 = np + ! first write to isomers.xyz and fragmentpairs.xyz and filter out + write (*,*) "writing isomers to and fragmentpairs to " + if (lprint) write (*,*) "writing product structures and molecular masses to directories" + + prstruc = env%msnfrag + if (prstruc .gt. 0) then + write (*,*) "Printing only ",prstruc," selected structures to products.xyz" + incr = np/prstruc + else + incr = 1 + prstruc = np + end if + + write (*,'(a)') '========================================================' + write (*,'(a,i0)') " number of printed structures: ",prstruc ! todo changes depending on print settings + write (*,*) " directory | fragment type | rel. energy [kcal/mol]" + if (.not.env%msiso) write (*,*) " fragmentpair: | sumformula | molecular mass" + + np2 = 0 ! number of printed sturctures + open (newunit=ich1,file=fname,status='replace') + do i = 1,np,incr + if (np2 .ge. prstruc) exit + erel = (mso%pl%mol(i)%etot-estart)*autokcal + if (mso%pl%mol(i)%fragcount .lt. 1) write (*,*) "no fragment count for ",i + if (mso%pl%mol(i)%fragcount .eq. 1) then + np2 = np2+1 + write (isodir,'(a,i0)') "p",np2 + write (*,'(a4,a20,f18.8)') trim(isodir),"isomer",erel + strucname = 'isomer.xyz' + elseif (mso%pl%mol(i)%fragcount .gt. 1) then + np2 = np2+1 + write (isodir,'(a,i0)') "p",np2 + write (*,'(a4,a20,f18.8)') trim(isodir),"fragmentpair",erel + strucname = 'pair.xyz' + end if + write (comment,'(f9.5,1x,a)') mso%pl%mol(i)%etot,trim(isodir) + call wrxyz(ich1,mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),trim(comment)) + + ! write to directories + if (lprint) then + r = makedir(trim(isodir)) + call chdir(trim(isodir)) + open (newunit=ich,file=strucname,status='replace') + call wrxyz(ich,mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),mso%pl%mol(i)%etot) + close (ich) + ! just write mass for pairs too + mass = 0 + do k = 1,nat + mass = mass+ams(mso%pl%mol(i)%at(k)) + end do + call wrshort_real("molmass",mass) + call chdir(trim(thisdir)) + end if + ! no write separated fragment structures + if (mso%pl%mol(i)%fragcount .gt. 1) then + do j = 1,2 ! loop over fragment 1 and 2 + natf(j) = count(mso%pl%mol(i)%fragi == j) + mass = 0 + atf(j,:) = 0 + do k = 1,nat + if (mso%pl%mol(i)%fragi(k) == j) then + mass = mass+ams(mso%pl%mol(i)%at(k)) + atf(j,k) = mso%pl%mol(i)%at(k) + end if + end do + sumformula = sumform(nat,atf(j,:)) + ! write to directories + if (lprint) then + write (fragdir,'(a,i0)') trim(isodir)//"f",j + r = makedir(trim(fragdir)) + call chdir(fragdir) + strucname = 'fragment.xyz' + open (newunit=ich,file=strucname,status='replace') + write (ich,*) natf(j) + write (ich,*) + do k = 1,nat + if (mso%pl%mol(i)%fragi(k) == j) then !only write xyz if fragment really exists + write (ich,'(a2,5x,3F18.8)') i2e(mso%pl%mol(i)%at(k)),mso%pl%mol(i)%xyz(1,k),mso%pl%mol(i)%xyz(2,k),mso%pl%mol(i)%xyz(3,k) + end if + end do + call wrshort_real("mass",mass) + call chdir(trim(thisdir)) + end if + !write(*,*) " directory | sumformula | mass" + write (*,'(a6,i0,5x,a20,9x,f9.5)') trim(isodir)//"f",j,trim(sumformula),mass + end do + end if + end do + + write (*,'(a)') '========================================================' + write (*,*) + call wrshort_int('npairs',np2) + if (env%msnoiso) then + write (*,*) "sorted out ",np0-np2,"non-dissociated structures" + elseif (env%msiso) then + write (*,*) "sorted out ",np0-np2,"dissociated structure pairs" + end if + write (*,*) "Number of generated isomers: ",nisomers + write (*,*) "Number of generated fragmentpairs: ",nfragpairs + end subroutine write_fragments + +!=======================================================================================! + + subroutine msinputreader(mso,msinput) +!********************************************************** +!* a simple input reader to read in additional input file +!* to set more advanced options +!********************************************************** + implicit none + type(msobj) :: mso + integer :: ich,io,iocheck + character(len=40) :: line + character(len=80) :: msinput + real(wp) :: xx(10) + integer :: j + logical :: ex + + write (*,*) "msinput is:,",trim(msinput) + !---- read input file + inquire (file=trim(msinput),exist=ex) + if (.not.ex) then + return + else + write (*,*) 'Reading <'//trim(msinput)//'> file' + end if + + mso%rcut = 1.3_wp !> cutoff for factor of sum of covalent radii for detection of fragmentation of fragment_structures + mso%T = 5000.0_wp !> electronic temperature! favours open-shell fragment generation + mso%fc = 0.05_wp !> repulsive force constant in constrained optimizations + mso%fc_attr = -0.05_wp !> attractive force constant in constrained optimizations + mso%cdist = 1.5_wp !> constraing distance scaling factor of rcov + mso%cdist_att = 0.5_wp !> constraing distance scaling factor of rcov for attractive part + mso%distthr_att = 4.0_wp !> distance threshold in Angstrom for attractive part of constraint + mso%atomshift = 0.75_wp !> shift of atoms in random atom displacement + mso%fragdist = 0.0_wp !> increase distance (in Angstrom) to separate fragments (useful for subsequent transition state search) + mso%maxc = 15 + + open (newunit=ich,file=trim(msinput),status='old') + + ! Read line-by-line + do + read (ich,'(a)',iostat=iocheck) line + + ! Check for errors in the input + if (iocheck > 0) then !Fail + write (*,*) 'Something is wrong in the input. Exiting...' + stop + + ! End-of-file + elseif (iocheck < 0) then !EOF + exit + + ! Keywords + else + write (*,'(''>'',a)') line + ! capatilize and trimthe line + line = uppercase(line) + line = trim(line) + + ! read keywords + if (index(line,'FRAGDIST') /= 0) then + call readl(trim(line),xx,j) + mso%fragdist = xx(1) + end if + if (index(line,'ATOMSHIFT') /= 0) then + call readl(trim(line),xx,j) + mso%atomshift = xx(1) + end if + if (index(line,'DISTTHR_ATTR') /= 0) then + call readl(trim(line),xx,j) + mso%distthr_att = xx(1) + end if + if (index(line,'FC_REP') /= 0) then + call readl(trim(line),xx,j) + mso%fc = xx(1) + end if + if (index(line,'FC_ATTR') /= 0) then + call readl(trim(line),xx,j) + mso%fc_attr = xx(1) + end if + if (index(line,'ETEMP') /= 0) then + call readl(trim(line),xx,j) + mso%T = xx(1) + end if + end if + + end do + + end subroutine msinputreader + +!=======================================================================================! + + subroutine detect_fragments(mso,lprint,nisomers,nfragpairs) +!******************************************************************* +! detection of fragments and removal of multiple fragmented species +!******************************************************************* + implicit none + type(msobj) :: mso + integer :: nat + integer :: i + integer :: nfragpairs,nisomers + integer :: npoly ! number of multiple fragmented species + integer :: np ! number of products + character(len=2048) :: multfrags + logical,intent(in) :: lprint + + write (*,*) "Detect fragments" + + nat = mso%pl%mol(1)%nat + np = mso%pl%nmol + + nisomers = 0 + nfragpairs = 0 + npoly = 0 + + ! determine fragments for every structure + do i = 1,np + allocate (mso%pl%mol(i)%fragi(nat)) + call fragment_structure(mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),mso%rcut,1,0,mso%pl%mol(i)%fragi,mso%pl%mol(i)%fragcount) ! works better than mrec + end do + ! remove multiple fragmented species + multfrags = '' + i = 0 + do + i = i+1 + if (i .gt. np) exit ! np changes within the loop + if (mso%pl%mol(i)%fragcount .gt. 2) then + npoly = npoly+1 + if (lprint) write (multfrags,'(a,i0)') trim(multfrags)//" ",i + call mso%pl%remove(i) + i = i-1 + np = np-1 + cycle + elseif (mso%pl%mol(i)%fragcount .eq. 2) then + nfragpairs = nfragpairs+1 + else ! for isomers do nothing + nisomers = nisomers+1 + cycle + end if + end do + if (lprint.and.len(trim(multfrags)) .gt. 0) then + write (*,'(a,i0,a)') "removed multiple fragmented structures: "//trim(multfrags) + end if + return + end subroutine detect_fragments + +!=======================================================================================! + + subroutine increase_fragdist(mso,lprint) +!**************************************************************************** +!* increase distance between fragments, better for transistion state search +!* we remove here also multiple fragmented species, +!* and unreasonable fragments (too much distance between fragments) +!**************************************************************************** + use axis_module + implicit none + type(msobj) :: mso + integer :: nat + integer,allocatable :: at1(:),at2(:) + integer,allocatable :: fragi(:) !fragment index + integer :: frag1,frag2,i,j + integer :: npoly ! number of multiple fragmented species + integer :: np ! number of products + real(wp) :: norm + real(wp),allocatable :: xyz1(:,:),xyz2(:,:) ! xyz of fragment 1 and 2 + real(wp) :: cmass1(3),cmass2(3) ! center of mass of fragment 1 and 2 + character(len=2048) :: distfrags + logical,intent(in) :: lprint + + frag1 = 0 + frag2 = 0 + if (mso%fragdist .gt. 0.0_wp) then + write (*,'(a,f10.8,a)') " Increase distance of fragments by ",mso%fragdist," Angstrom" + else + return + end if + + nat = mso%pl%mol(1)%nat + np = mso%pl%nmol + + ! increase distance between fragments (based on center of mass) and remove unreasonable fragments + ! that are too far away from each other, which cause problems in transition state search + + distfrags = '' + i = 0 + do + i = i+1 + if (i .gt. np) exit ! np changes within the loop + ! increase distance between fragments + if (mso%pl%mol(i)%fragcount .eq. 2) then + do j = 1,nat + if (mso%pl%mol(i)%fragi(j) == 1) then + frag1 = frag1+1 + end if + if (mso%pl%mol(i)%fragi(j) == 2) then + frag2 = frag2+1 + end if + end do + allocate (at1(frag1),at2(frag2)) + allocate (xyz1(3,frag1),xyz2(3,frag2)) + frag1 = 0 + frag2 = 0 + do j = 1,nat + if (mso%pl%mol(i)%fragi(j) == 1) then + frag1 = frag1+1 + xyz1(:,frag1) = mso%pl%mol(i)%xyz(:,j) + at1(frag1) = mso%pl%mol(i)%at(j) + end if + if (mso%pl%mol(i)%fragi(j) == 2) then + frag2 = frag2+1 + xyz2(:,frag2) = mso%pl%mol(i)%xyz(:,j) + at2(frag2) = mso%pl%mol(i)%at(j) + end if + end do + call CMAv(frag1,at1,xyz1,cmass1) + call CMAv(frag2,at2,xyz2,cmass2) + norm = sqrt((cmass2(1)-cmass1(1))**2+(cmass2(2)-cmass1(2))**2+(cmass2(3)-cmass1(3))**2) + ! remove structures with unreasonable distance between fragments + if (norm .gt. 15) then !15 angstrom worked well in tests + if (lprint) write (distfrags,'(a,i0)') trim(distfrags)//" ",i + call mso%pl%remove(i) + i = i-1 + np = np-1 + cycle + end if + do j = 1,nat + ! move fragment 2 away from fragment 1 + if (mso%pl%mol(i)%fragi(j) == 2) then + mso%pl%mol(i)%xyz(:,j) = mso%pl%mol(i)%xyz(:,j)+(cmass2-cmass1)/norm*mso%fragdist + end if + end do + deallocate (at1,at2,xyz1,xyz2) + else ! for isomers do nothing + cycle + end if + end do + + if (lprint.and.len(trim(distfrags)) .gt. 0) then + write (*,'(a,i0,a)') "removed unreasonable structures:"//trim(distfrags) + end if + + return + end subroutine increase_fragdist + +!=======================================================================================! + + subroutine rdplist(mso,estart,fname) +!******************************************************* +!* collect structures of ensemble fil into product list +!******************************************************* + use atmasses + implicit none + integer :: nat + integer :: nall + integer :: chrg ! charge of the molecule + type(msobj) :: mso + character(len=*) :: fname + integer,allocatable :: at(:) + real(wp),allocatable :: xyz(:,:) + real(wp) :: eread + real(wp) :: molmass + real(wp) :: estart + character(len=6) :: sym + character(len=512) :: line + integer :: i,j,k,dum,io,ich + logical :: ex,found + + found = .false. + write (*,*) "reading product structures file ",fname + + call rdensembleparam(fname,nat,nall) + allocate (at(nat),xyz(3,nat)) + call rdshort('.CHRG',chrg) + eread = 0.0_wp + xyz = 0.0_wp + open (newunit=ich,file=fname) + do i = 1,nall + read (ich,*,iostat=io) dum + if (io < 0) exit + if (io > 0) cycle + read (ich,'(a)',iostat=io) line + if (io < 0) exit + eread = grepenergy(line) + if (abs(eread-estart) .lt. 0.00001_wp) then ! skip input structure, detect by energy + found = .true. + do k = 1,dum + read (ich,'(a)',iostat=io) line + end do + cycle + if (io < 0) exit + end if + do j = 1,dum + read (ich,'(a)',iostat=io) line + if (io < 0) exit + call coordline(line,sym,xyz(1:3,j),io) + if (io .ne. 0) then + backspace (ich) + exit + end if + at(j) = e2i(sym) + end do + molmass = molweight(nat,at) + call mso%pl%append(nat,at,xyz,eread,chrg,0,molmass) + end do + close (ich) + + if (io < 0) then + error stop 'error while reading product file.' + end if + + deallocate (xyz,at) + if (.not.found) then + write (*,*) "Warning, input structure not found in ensemble and could not be sorted out" + end if + return + end subroutine rdplist + +!=======================================================================================! + + subroutine wrplist(mso,fname) +!************************* +!* write products to file +!************************* + implicit none + integer :: i,ich,io + character(len=*) :: fname + type(msobj) :: mso + open (newunit=ich,file=fname) + do i = 1,mso%pl%nmol + call wrxyz(ich,mso%pl%mol(i)%nat,mso%pl%mol(i)%at,mso%pl%mol(i)%xyz(:,:),mso%pl%mol(i)%etot) + end do + close (ich) + return + end subroutine wrplist + +!=======================================================================================! +!=======================================================================================! +end module msmod !> END OF MODULE +!=======================================================================================! +!=======================================================================================! diff --git a/src/msreact/msreact.f90 b/src/msreact/msreact.f90 new file mode 100644 index 00000000..c3b3e857 --- /dev/null +++ b/src/msreact/msreact.f90 @@ -0,0 +1,556 @@ +!================================================================================! +! This file is part of crest. +! +! Copyright (C) 2020 Stefan Grimme, Philipp Pracht +! +! crest 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 3 of the License, or +! (at your option) any later version. +! +! crest 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 crest. If not, see . +!================================================================================! + +!==============================================================! +! the handler function is the sub-program that is called +! by crest, i.e., within this routine the seperate input file is +! read and the calculation is started. +!==============================================================! +subroutine msreact_handler(env,tim) + use crest_parameters + use crest_data + use msmod + use strucrd + use atmasses + use iomod + use zdata + implicit none + + type(systemdata) :: env + type(timer) :: tim + + type(msobj) :: mso !a new msreact object + type(coord) :: struc + + interface + subroutine msreact_topowrap(mol,pair,paths,wboname) + import :: msmol + implicit none + type(msmol) :: mol + integer :: pair(mol%nat*(mol%nat+1)/2) + integer :: paths(mol%nat*(mol%nat+1)/2,mol%nat) + character(len=*),optional :: wboname + end subroutine msreact_topowrap + end interface + integer :: nat + integer,allocatable :: pair(:) + integer,allocatable :: paths(:,:) + integer :: k + integer :: nisomers,nfragpairs + real(wp) :: molmass + real(wp) :: estart ! energy of input structure + + call tim%start(1,'MSREACT') + + call msreact_head(env,mso) + + ! get starting energy of lowest fragment + call get_input_energy(env,mso%T,estart) + !-- read the input coord and put it into the + ! product-list as first structure + call struc%open('coord') + struc%xyz = struc%xyz*bohr !to Angstrom, from this point on by convention! + molmass = molweight(struc%nat,struc%at) + call mso%pl%append(struc%nat,struc%at,struc%xyz,estart,env%chrg,1,molmass) + if (env%mslargeprint) call wrshort_real("mass",mso%pl%mol(1)%molmass) ! write mass + call struc%deallocate() + + !-- additional input file can be read here + call msinputreader(mso,env%msinput) + + nat = mso%pl%mol(1)%nat + k = nat*(nat+1)/2 + allocate (pair(k),paths(k,nat)) + call msreact_topowrap(mso%pl%mol(1),pair,paths,'wbo') + + !-- setting the threads for correct parallelization + if (env%autothreads) then + call ompautoset(env%threads,6,env%omp,env%MAXRUN,0) !set global OMP/MKL variable for xtb jobs + end if + + !-- do the directory setup and optimizations + call msreact(env,mso,nat,pair) + deallocate (paths,pair) + + if (env%mslargeprint) then + write (*,*) "All structures written to file " + call wrplist(mso,'crest_allproducts.xyz') + end if + ! detect fragmented structures + call detect_fragments(mso,env%mslargeprint,nisomers,nfragpairs) + call increase_fragdist(mso,env%mslargeprint) + write (*,*) "Remaining structures:",mso%pl%nmol-1 + + ! call write_products(env,mso) print out structures if lprint is chosen? + + ! sort out topologically identical structures + call sortoutduplicates(env,mso) + + write (*,'(a,f10.4,a)') "Calling CREGEN routine and sort out structures above energy window of",env%ewin," kcal/mol" + ! write plist to file to use cregen routine + call wrplist(mso,'crest_unique_products.xyz') + if (env%mslargeprint) then + write (*,*) "All unique structures written to file " + end if + ! sort according to energies and RMSD, remove highe energy structures + call newcregen(env,13,'unique_products.xyz') + call mso%pl%dealloc() + ! read in again plist after cregen and remove starting structure + call rdplist(mso,estart,'crest_msreact_products.xyz') + write (*,*) "Remaining structures after cregen:",mso%pl%nmol + call detect_fragments(mso,env%mslargeprint,nisomers,nfragpairs) + + ! write final structures to file and directories + call write_fragments(env,mso,estart,nisomers,nfragpairs,'crest_msreact_products.xyz',env%mslargeprint) + + !cleanup files + if (.not.env%mslargeprint) then + call rmrf('MSDIR') + call rmrf('crest_allproducts.xyz') + call rmrf('crest_unique_products.xyz') + end if + call rmrf('crest_unique_products.sorted') + call rmrf('scoord.1') + call rmrf('xtbscreen.xyz xtblmoinfo') + call tim%stop(1) + return +end subroutine msreact_handler + +!==============================================================! +! the main implementation of the msreact algo +!==============================================================! +subroutine msreact(env,mso,nat,pair) + use crest_parameters + use msmod + use iomod + use miscdata,only:rcov + use utilities,only:lin + implicit none + + type(msobj) :: mso !main storage object + type(msmol) :: inmol ! xyz of molecule + type(msmol) :: inmoldisp ! xyz of molecule with displacement + type(systemdata) :: env ! system data + + integer :: nat + integer :: pair(nat*(nat+1)/2) ! gives number of bonds between two atoms lin(i,j) + integer :: nbonds + integer :: i,j,k,l + integer :: p + integer :: np + integer :: io + integer :: nbaseat ! number of lewis basic atoms according to a xTB LMO analysis + integer,allocatable :: basicatlist(:) ! list of lewis basic atoms according to a xTB LMO analysis + + character(len=:),allocatable :: subdir + character(len=40) :: pdir + character(len=30) :: base + character(len=512) :: thisdir + real(wp) :: constr_dist + + inmol = mso%pl%mol(1) + + !-- main subdirectory handling + call getcwd(thisdir) + subdir = 'MSDIR' + io = makedir(subdir) + call chdir(subdir) + + !-- get specific pairs + base = 'Pair_' + np = 0 + do i = 1,nat + do j = i,nat + k = lin(i,j) + p = pair(k) + if (p .eq. 1.and.(inmol%at(i) .eq. 1.or.inmol%at(j) .eq. 1)) cycle ! do not distort all X-H bonds + if (p .gt. 0.and.p .le. env%msnbonds) then + !write(*,*) " distort bond ",i,"-",j," with ",p," bonds" + np = np+1 + write (pdir,'(a,i0)') trim(base),np + constr_dist = mso%cdist*(rcov(inmol%at(i))+rcov(inmol%at(j)))*bohr+float(p) + ! write(*,*) inmol%at(i),inmol%at(j),constr_dist + call isodir(mso,trim(pdir),inmol,i,j,constr_dist,mso%fc) + end if + end do + end do + + !attractive potential for H-shifts + !for H-shifts: xtb -lmo , read xtb.out -> get list of atoms with pi or lone-pair -> H-allowed to shift there + if (env%msattrh) then + call chdir(thisdir) + call xtblmo(env,.true.) + call readbasicpos(env,nbaseat,basicatlist) + call chdir(subdir) + write (*,*) "Add attractive potential for H-shifts:" + !-- setting the threads for correct parallelization (were changed in xtblmo call) + if (env%autothreads) then + call ompautoset(env%threads,6,env%omp,env%MAXRUN,0) !set global OMP/MKL variable for xtb jobs + end if + do i = 1,nat + do j = i,nat + k = lin(i,j) + p = pair(k) + if (p .ge. 2) then ! only add attractive potential for pairs with at least 2 bonds distance + if ((inmol%at(i) .eq. 1.and.(findloc(basicatlist,j,1) /= 0)).or. & ! pair of H and basic atom? + & (inmol%at(j) .eq. 1.and.(findloc(basicatlist,i,1) /= 0))) then + if (msmoldist(inmol,i,j) .lt. mso%distthr_att) then + np = np+1 + !write(*,*) "add attractive potential for pair ", np," between ",i," and ",j + write (pdir,'(a,i0)') 'Pair_',np + constr_dist = mso%cdist_att*(rcov(inmol%at(i))+rcov(inmol%at(j)))*bohr+float(p) + call isodir(mso,trim(pdir),inmol,i,j,constr_dist,mso%fc_attr) + end if + end if + end if + end do + end do + end if + + write (*,*) '# of distortions',np + call msreact_jobber(env,np,base,.true.,.false.) + call msreact_collect(nat,np,base,mso) + + !-- do additional shifting of atoms and subsequent optimization + if (env%msnshifts .gt. 0) then + write (*,'(a,i0,a)') "Shift atoms ",env%msnshifts," times and reoptimize" + base = 'Distortion_' + np = 0 + do l = 1,env%msnshifts + call shiftatoms(mso,inmol,inmoldisp) + np = np+1 + write (pdir,'(a,i0)') trim(base),np + call isodiropt(mso,trim(pdir),inmoldisp) + end do + call msreact_jobber(env,np,base,.false.,.false.) + ! collect everything in isomer list + call msreact_collect(nat,np,base,mso) + end if + + ! do additional shifting of atoms and subsequent optimization with constrained bonds + if (env%msnshifts2 .gt. 2) then + write (*,'(a,i0,a)') "Shift atoms",env%msnshifts2," times and apply repulsive potential on bonds of distorted structure" + base = 'Distortedpair_' + np = 0 + do l = 1,env%msnshifts2 + call shiftatoms(mso,inmol,inmoldisp) + do i = 1,nat + do j = i,nat + k = lin(i,j) + p = pair(k) + if (p .eq. 1.and.(inmol%at(i) .eq. 1.or.inmol%at(j) .eq. 1)) cycle + if (pair(k) .le. env%msnbonds) then + np = np+1 + write (pdir,'(a,i0)') trim(base),np + constr_dist = mso%cdist*(rcov(inmoldisp%at(i))+rcov(inmoldisp%at(j)))*bohr+float(p) + call isodir(mso,trim(pdir),inmoldisp,i,j,constr_dist,mso%fc) + end if + end do + end do + end do + call msreact_jobber(env,np,base,.true.,.false.) + ! collect everything in isomer list + call msreact_collect(nat,np,base,mso) + end if + + call chdir(thisdir) + return +end subroutine msreact + +!============================================================! +! make a dir for a constrained optimization of a structure, +! a controlfile with constraints on atoms A and B (at dist D) +! will be written into the directory +!============================================================! +subroutine isodir(mso,dirname,mol,A,B,D,fc) + use crest_parameters + use msmod + use iomod + use strucrd,only:wrxyz + implicit none + type(msobj) :: mso + character(len=*) :: dirname + type(msmol),intent(in) :: mol + integer :: A,B + real(wp) :: D + real(wp),intent(in) :: fc ! force constant for bias potential + + character(len=:),allocatable :: fname + character(len=20) :: dumm + integer :: io,ich + + io = makedir(dirname) !create the directory + + fname = trim(dirname)//'/'//'struc.xyz' + open (newunit=ich,file=fname) + call wrxyz(ich,mol%nat,mol%at,mol%xyz) + close (ich) + + fname = trim(dirname)//'/'//'.CHRG' + open (newunit=ich,file=fname) + write (ich,'(i0)') mol%chrg ! EI +1, DEA -1, CID 0 + close (ich) + + ! write control file for bond length constraint + + fname = trim(dirname)//'/'//'.xc1' + open (newunit=ich,file=fname) + write (ich,'(a)') '$scc' + write (dumm,'(f16.2)') mso%T + write (ich,'(1x,a,a)') 'temp=',adjustl(trim(dumm)) + write (ich,'(a)') '$constrain' + write (dumm,'(f16.4)') fc + write (ich,'(3x,a,a)') 'force constant=',adjustl(trim(dumm)) + write (ich,'(3x,a,1x,i0,a,1x,i0,a,1x,f8.5)') 'distance:',A,',',B,',',D + close (ich) + + fname = trim(dirname)//'/'//'.xc2' + open (newunit=ich,file=fname) + write (ich,'(a)') '$scc' + write (dumm,'(f16.2)') mso%T + write (ich,'(1x,a,a)') 'temp=',adjustl(trim(dumm)) + write (ich,'(a)') '$opt' + write (ich,'(1x,a)') 'maxcycle=',mso%maxc ! only allow 15 cycles to avoid recombination of fragments + write (ich,'(a)') '$write' + write (ich,'(1x,a)') 'wiberg=true' + close (ich) + + return +end subroutine isodir + +!============================================================! +! make a dir for a optimization of a displaced structure, +! a controlfile for an optimization after the atom shifting +! will be written into the directory +!============================================================! +subroutine isodiropt(mso,dirname,mol) + use iso_fortran_env,only:wp => real64 + use msmod + use iomod + use strucrd,only:wrxyz + implicit none + type(msobj) :: mso + character(len=*) :: dirname + type(msmol),intent(in) :: mol + integer :: A,B + real(wp) :: D + + character(len=:),allocatable :: fname + character(len=20) :: dumm + integer :: io,ich + + io = makedir(dirname) !create the directory + + fname = trim(dirname)//'/'//'struc.xyz' + open (newunit=ich,file=fname) + call wrxyz(ich,mol%nat,mol%at,mol%xyz) + close (ich) + + fname = trim(dirname)//'/'//'.CHRG' + open (newunit=ich,file=fname) + write (ich,'(i0)') mol%chrg ! EI +1, DEA -1, CID 0, give in crest call + close (ich) + + fname = trim(dirname)//'/'//'.xc1' + open (newunit=ich,file=fname) + write (ich,'(a)') '$scc' + write (dumm,'(f16.2)') mso%T + write (ich,'(1x,a,a)') 'temp=',adjustl(trim(dumm)) + close (ich) + return +end subroutine isodiropt + +!=====================================================================! +! The job construction routine for MSREACT +! (will have to be modified later, for now it is for testing) +!=====================================================================! +subroutine msreact_jobber(env,ndirs,base,constr,niceprint) + use crest_parameters + use msmod + use iomod + implicit none + integer :: ndirs + character(len=*) :: base + logical :: niceprint,constr + character(len=1024) :: jobcall,jobcall2 + type(systemdata) :: env ! system data + + jobcall = '' + jobcall2 = '' + + if (constr) then + write (jobcall,'(a)') 'xtb struc.xyz --opt loose --input .xc1 '//trim(env%gfnver)//' > split.out 2>/dev/null' + write (jobcall2,'(a)') 'xtb xtbopt.xyz --opt crude --input .xc2 '//trim(env%gfnver)//' > xtb.out 2>/dev/null' + jobcall = trim(jobcall)//' ; '//trim(jobcall2) + else + write (jobcall,'(a)') 'xtb struc.xyz --opt crude --input .xc1 '//trim(env%gfnver)//' > xtb.out 2>/dev/null' + end if + + !-- directories must be numbered consecutively + call opt_OMP_loop(ndirs,base,jobcall,niceprint) + write (*,*) + write (*,*) 'done.' + return +end subroutine msreact_jobber + +!=====================================================================! +! A wrapper to generate the topology for a molecule within the +! MSREACT subprogram +!=====================================================================! +subroutine msreact_topowrap(mol,pair,paths,wboname) + use crest_parameters + use msmod + use zdata + use adjacency + use utilities,only:lin + implicit none + type(msmol) :: mol + integer :: pair(mol%nat*(mol%nat+1)/2) + !integer :: pair(mol%nat,mol%nat) + integer :: paths(mol%nat*(mol%nat+1)/2,mol%nat) + character(len=*),optional :: wboname + type(zmolecule) :: zmol + + integer,allocatable :: A(:,:) + integer,allocatable :: prev(:,:) + real(wp),allocatable :: E(:,:) + real(wp),allocatable :: dist(:,:) + + integer :: lpath,i,j,k + integer,allocatable :: path(:) + logical :: ex + + ex = .false. + if (present(wboname)) then + inquire (file=wboname,exist=ex) + end if + if (ex) then + call simpletopo(mol%nat,mol%at,mol%xyz,zmol,.false.,.false.,wboname) + else + mol%xyz = mol%xyz/bohr !CN based topo requires Bohrs + call simpletopo(mol%nat,mol%at,mol%xyz,zmol,.false.,.false.,'') + mol%xyz = mol%xyz*bohr + end if + + allocate (A(mol%nat,mol%nat),E(mol%nat,mol%nat)) + call zmol%adjacency(A,E) + + allocate (prev(mol%nat,mol%nat),dist(mol%nat,mol%nat)) + + call FloydWarshall(mol%nat,A,E,dist,prev) + allocate (path(mol%nat),source=0) + do i = 1,mol%nat + do j = i,mol%nat + path = 0 + call getPathFW(mol%nat,prev,i,j,path,lpath) + !write(*,*) path(1:lpath) + k = lin(i,j) + pair(k) = lpath-1 ! number of bonds + paths(k,:) = path(:) + end do + end do + + deallocate (dist,prev) + deallocate (E,A) + + call zmol%deallocate() !clear the zmol memory + return +end subroutine msreact_topowrap + +!========================================================================! +! collect structures of optimized molecules +! xyz files should still have the same number and order of atoms +!========================================================================! +subroutine msreact_collect(nat,np,base,mso) + use crest_parameters + use strucrd + use msmod + use atmasses + implicit none + integer :: nat + integer :: np + integer :: ich + integer :: chrg ! charge of the molecule + type(msobj) :: mso + character(len=40) :: pdir + character(len=*) :: base + character(len=:),allocatable :: optfile + character(len=128) :: newcomment + integer :: p + logical :: ex + integer,allocatable :: at(:) + real(wp),allocatable :: xyz(:,:) + real(wp) :: etot + real(wp) :: molmass + + allocate (at(nat),xyz(3,nat)) + + do p = 1,np + write (pdir,'(a,i0)') trim(base),p + optfile = trim(pdir)//'/'//'xtbopt.xyz' + inquire (file=optfile,exist=ex) + if (ex) then + call rdcoord(optfile,nat,at,xyz,etot) + xyz = xyz*bohr + call rdshort(trim(pdir)//'/'//'.CHRG',chrg) + molmass = molweight(nat,at) + call mso%pl%append(nat,at,xyz,etot,chrg,0,molmass) + end if + end do + close (ich) + + deallocate (xyz,at) + return +end subroutine msreact_collect + +!============================================================! +! shift atoms of input structure randomly +! to generate more structures +! for planar molecules important +!============================================================! +subroutine shiftatoms(mso,molin,molout) + use iso_fortran_env,only:wp => real64 + use msmod + use iomod + use strucrd,only:wrxyz + implicit none + type(msobj) :: mso + type(msmol) :: molin + type(msmol) :: molout + real(wp) :: D + real(wp) :: x,shift + real(wp) :: dist ! distance to shift atoms + character(len=80) :: fname + character(len=20) :: dumm + integer :: io,ich,i,j,count,incr + + dist = mso%atomshift ! distance of displacement + molout = molin + do i = 1,molin%nat + do j = 1,3 + call Random_Number(x) + ! positive and negative direction possible + x = 2.0_wp*x-1.0_wp + shift = x*dist + molout%xyz(j,i) = molin%xyz(j,i)+shift + end do + end do + return +end subroutine shiftatoms + diff --git a/src/printouts.f90 b/src/printouts.f90 index bd124fb7..7d738c20 100644 --- a/src/printouts.f90 +++ b/src/printouts.f90 @@ -299,7 +299,34 @@ subroutine confscript_morehelp(flag) write (*,'(5x,''-freqscal : defines frequency scale factor. Only for outprint'')') write (*,'(5x,''-freqlvl [method] : define a method for frequency computation. All gfn versions are supported'')') write (*,*) - + + case('msreact') + write (*,'(1x,'' mass spectral fragment generator (msreact)'')') + write (*,'(1x,''General usage :'')') + write (*,'(5x,'' -msreact [options]'')') + write (*,'(1x,''options:'')') + write(*,'(5x,''-msnoattrh : deactivate attractive potential between hydrogen and LMO centers)'')') + write(*,'(5x,''-msnshifts [int] : perform n optimizations with randomly shifted atom postions (default 0) '')') + write(*,'(5x,''-msnshifts2 [int] : perform n optimizations with randomly shifted atom postions and repulsive potential applied to bonds (default 0) '')') + write(*,'(5x ''-msnbonds [int] : maximum number of bonds between atoms pairs for applying repulsive potential (default 3)'')') + write(*,'(5x,''-msmolbar : sort out topological duplicates by molbar codes (requires sourced "molbar")'')') + write(*,'(5x,''-msinchi : sort out topological duplicates by inchi codes (requires sourced "obabel")'')') + write(*,'(5x ''-msnfrag [int] : number of fragments that are printed by msreact (random selection)'')') + write(*,'(5x,''-msiso : print only non-dissociated structures (isomers)'')') + write(*,'(5x,''-msnoiso : print only dissociated structures'')') + write(*,'(5x,''-mslargeprint : do not remove temporary files and MSDIR do not remove temporary files and MSDIR with constrained optimizations'')') + write (*,'(5x,''-chrg : set the molecules´ charge'')') + write (*,'(5x,''-ewin : set energy window in for sorting out fragments kcal/mol,'')') + write (*,'(5x,'' [default: 200.0 kcal/mol] '')') + write (*,'(5x,''-msinput : read in an input file with special settings for msreact'')') + write (*,'(5x,''keywords for inputfile:'')') + write (*,'(5x,'' fragdist : increase distance between fragments xyz structures (default 0 Angstrom) '')') + write (*,'(5x,'' atomshift : shift of atoms in random atom displacement (default 0.75 Angstrom)'')') + write (*,'(5x,'' distthr_attr : distance threshold in Angstrom for H-LMO attraction (default 4.0 Angstrom) '')') + write (*,'(5x,'' fc_rep : force constant for repulsive potential between atom pairs (default 0.5) '')') + write (*,'(5x,'' fc_attr : force constant for attractive potential between hydrogen and LMO centers (default -0.5) '')') + write (*,'(5x,'' etemp : electronic temperature in xTB optimizations'')') + case ('other') write (*,'(1x,''Other tools for standalone use:'')') write (*,'(5x,''-zsort : use only the zsort subroutine'')') @@ -504,6 +531,27 @@ end subroutine qcg_head !========================================================================================! +subroutine msreact_head() + implicit none + write (*,*) + write (*,'(2x,''========================================'')') + write (*,'(2x,''| |'')') + write (*,'(2x,''| MSREACT |'')') + write (*,'(2x,''| automated MS fragment generator |'')') + write (*,'(2x,''| |'')') + write (*,'(2x,''| University of Bonn, MCTC |'')') + write (*,'(2x,''========================================'')') + write (*,'(2x,'' S. Grimme, P. Pracht, J. Gorges.'')') + write (*,*) + write (*,'(3x,''Cite work conducted with this code as'')') + write (*,'(/,3x,''Philipp Pracht, Stefan Grimme, Christoph Bannwarth, Fabian Bohle, Sebastian Ehlert, Gereon Feldmann,'')') + write (*,'(3x,''Johannes Gorges, Marcel Müller, Tim Neudecker, Christoph Plett, Sebastian Spicher, Pit Steinbach,'')') + write (*,'(3x,''Patryk A. Wesolowski, and Felix Zeller J. Chem. Phys., 2024, submitted.'')') + write (*,*) + end subroutine msreact_head + +!========================================================================================! + subroutine smallhead(str) !********************************************** !> convert a string in a small header printout @@ -671,7 +719,7 @@ subroutine print_crest_metadata() write (*,'(2x,a,1x,a)') '-DWITH_GFNFF :',gfnffvar write (*,'(2x,a,1x,a)') '-DWITH_TBLITE :',tblitevar write (*,'(2x,a,1x,a)') '-DWITH_XHCFF :',xhcffvar - + write (*,'(2x,a,1x,a)') '-DWITH_LWONIOM :',lwoniomvar end subroutine print_crest_metadata !========================================================================================!