Merge pull request #483 from brunostega/main

added resdata and CM pbc warning

tools/cmdata/src/cmdata/cmdata.hpp

@@ -375,7 +375,11 @@ class CMData
     }
     if(!check_same && same_) same_ = false;
     if(natmol2_.size()<2) cross_ = false; 
-
+    if(nindex_>1 && (same_ || cross_))
+    {
+      printf("\n\n::::::::::::WARNING::::::::::::\nMore than 1 molcule found in the system.\nFix pbc before running cmdata using pbc mol\n");
+      printf(":::::::::::::::::::::::::::::::\n\n");
+    }
     if (same_)
     {
       f_inter_mol_same_ = cmdata::density::inter_mol_same_routine;

tools/resdata/.gitignore

@@ -0,0 +1,2 @@
+build/
+exe/

tools/resdata/CMakeLists.txt

@@ -0,0 +1,62 @@
+cmake_minimum_required(VERSION 3.16)
+
+project(resdata VERSION 0.1
+               DESCRIPTION "A programm to calculate contact data from gromacs trajectories for multi-eGO"
+               LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 17)
+set(RESDATA resdata)
+
+set(CMAKE_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
+set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
+
+add_executable(${RESDATA} main.cpp)
+
+# set build type to release if not set (check if its debug, Debug, or DEBUG)
+if(CMAKE_BUILD_TYPE MATCHES "Debug")
+  set(CMAKE_BUILD_TYPE Debug)
+  set(CMAKE_CXX_FLAGS_DEBUG "-g")
+elseif(CMAKE_BUILD_TYPE MATCHES "Release" OR NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+  set(CMAKE_CXX_FLAGS "-Wall -Wextra -march=native -Wno-unused-parameter")
+  set(CMAKE_CXX_FLAGS_RELEASE "-O3")
+endif()
+
+# find gromacs
+find_package(GROMACS REQUIRED NAMES gromacs gromacs_mpi gromacs_d gromacs_mpi_d HINTS "$ENV{GROMACS_DIR}")
+
+# include source and header files
+target_include_directories(${RESDATA} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/src)
+
+include(FetchContent)
+SET(FETCHCONTENT_QUIET OFF)
+SET(FETCHCONTENT_BASE_DIR ${CMAKE_CURRENT_BINARY_DIR}/resdata_fetch)
+
+message(STATUS "Fetching popt")
+FetchContent_Declare(
+  popt
+  URL https://github.com/rpm-software-management/popt/archive/refs/heads/master.zip
+  GIT_TAG 2bca0aa
+  PATCH_COMMAND patch --directory=${FETCHCONTENT_BASE_DIR}/popt-src -p0 < ${CMAKE_CURRENT_SOURCE_DIR}/popt.patch
+)
+FetchContent_MakeAvailable(popt)
+
+message(STATUS "Fetching xdrfile")
+FetchContent_Declare(
+  xdrfile
+  URL https://github.com/multi-ego/xdrfile/archive/refs/heads/chemfiles.zip 
+)
+FetchContent_MakeAvailable(xdrfile)
+
+# link libraries
+target_link_libraries(${RESDATA} PRIVATE Gromacs::libgromacs xdrfile popt)
+
+set_target_properties(${PROJECT_NAME} PROPERTIES INSTALL_RPATH ${CMAKE_INSTALL_PREFIX}/lib)
+install(TARGETS ${RESDATA} DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
+
+# build test
+if(RESDATA_BUILD_TESTS)
+  enable_testing()
+  add_subdirectory(test/)
+endif()

tools/resdata/README.md

@@ -0,0 +1,36 @@
+# resdata
+
+`resdata` is a tool that calculates residue wise interatomic distances and contact probabilities probabilities accounting for symmetric molecules . It works using the GROMACS API.
+
+## Installation
+
+To install it you need to compile and install GROMACS (2023/2024) using `-DGMX_INSTALL_LEGACY_API=ON` `-DBUILD_SHARED_LIBS=ON` 
+`-DGMX_INSTALL_NBLIB_API=ON`. Then in a local subfolder (e.g., `build`, run `cmake ..` then `make` and `make install`).
+
+## Usage
+resdata calculates a residue wise contact and distance matrix for complex system.
+It can Calculate block averages of probabilities
+It takes in input an index (only hole molecules can be removed otherwise an error is raised)
+
+```
+Usage: cmdata [OPTION...]
+  -f, --traj=FILE               Trajectory file
+  -s, --top=FILE                Topology file
+  -b, --t_begin[=FLOAT]         Start time
+  -e, --t_end[=FLOAT]           End time
+  -o, --out[=STRING]            Output prefix
+  -n, --index[=STRING]          Index file
+      --dt[=INT]                Time step
+      --cutoff[=DOUBLE]         Cutoff distance
+      --mol_cutoff[=DOUBLE]     Molecule cutoff distance
+      --nskip[=INT]             Number of frames to skip
+      --num_threads[=INT]       Number of threads
+      --mode[=STRING]           Mode of operation
+      --weights[=FILE]          Weights file
+      --no_pbc                  Ignore pbcs
+      --blk_avg                 Calculates block average for probabili
+
+Help options:
+  -?, --help                    Show this help message
+      --usage                   Display brief usage message
+```

tools/resdata/main.cpp

@@ -0,0 +1,169 @@
+// standard library imports
+#include <iostream>
+#include <string>
+#include <filesystem>
+// cmdata import
+#include "src/resdata/resdata.hpp"
+// external library import
+#include <popt.h>
+
+int main(int argc, const char** argv)
+{
+  std::cout << "################################################" << std::endl;
+  std::cout << "#################### RESDATA ####################" << std::endl;
+  std::cout << "################################################" << std::endl;
+  std::cout << "################## Version 0 ################" << std::endl;
+  std::cout << "################################################\n" << std::endl;
+
+  double cutoff = 0.55, mol_cutoff = 6.0;
+  int nskip = 0, num_threads = 1, mol_threads = -1, dt = 0;
+  float t_begin = 0.0, t_end = -1.0;
+  char *p_traj_path = NULL, *p_top_path = NULL, *p_mode = NULL, *p_weights_path = NULL, *p_index_path = NULL;
+  char *p_out_prefix = NULL;
+  std::string traj_path, top_path, mode, weights_path, index_path;
+  std::string out_prefix;
+  int *p_nopbc = NULL;
+  int *p_blk_avg = NULL;
+  bool nopbc = false;
+  bool blk_avg = false;
+
+  // make popt options
+  struct poptOption optionsTable[] = {
+    POPT_AUTOHELP
+    {"traj",        'f',  POPT_ARG_STRING,                          &p_traj_path,     0, "Trajectory file",             "FILE"},
+    {"top",         's',  POPT_ARG_STRING,                          &p_top_path,      0, "Topology file",               "FILE"},
+    {"t_begin",     'b',  POPT_ARG_FLOAT | POPT_ARGFLAG_OPTIONAL,   &t_begin,         0, "Start time",                  "FLOAT"},
+    {"t_end",       'e',  POPT_ARG_FLOAT | POPT_ARGFLAG_OPTIONAL,   &t_end,           0, "End time",                    "FLOAT"},
+    {"out",         'o',  POPT_ARG_STRING | POPT_ARGFLAG_OPTIONAL,  &p_out_prefix,    0, "Output prefix",               "STRING"},
+    {"index",       'n',  POPT_ARG_STRING | POPT_ARGFLAG_OPTIONAL,  &p_index_path,    0, "Index file   ",               "FILE"},
+    {"dt",          '\0', POPT_ARG_INT | POPT_ARGFLAG_OPTIONAL,     &dt,              0, "Time step",                   "INT"},
+    {"cutoff",      '\0', POPT_ARG_DOUBLE | POPT_ARGFLAG_OPTIONAL,  &cutoff,          0, "Cutoff distance",             "DOUBLE"},
+    {"mol_cutoff",  '\0', POPT_ARG_DOUBLE | POPT_ARGFLAG_OPTIONAL,  &mol_cutoff,      0, "Molecule cutoff distance",    "DOUBLE"},
+    {"nskip",       '\0', POPT_ARG_INT | POPT_ARGFLAG_OPTIONAL,     &nskip,           0, "Number of frames to skip",    "INT"},
+    {"num_threads", '\0', POPT_ARG_INT | POPT_ARGFLAG_OPTIONAL,     &num_threads,     0, "Number of threads",           "INT"},
+    {"mol_threads", '\0', POPT_ARG_INT | POPT_ARGFLAG_OPTIONAL,     &mol_threads,     0, "Number of molecule threads",  "INT"},
+    {"mode",        '\0', POPT_ARG_STRING | POPT_ARGFLAG_OPTIONAL,  &p_mode,          0, "Mode of operation",           "STRING"},
+    {"weights",     '\0', POPT_ARG_STRING | POPT_ARGFLAG_OPTIONAL,  &p_weights_path,  0, "Weights file",                "FILE"},
+    {"no_pbc",      '\0', POPT_ARG_NONE | POPT_ARGFLAG_OPTIONAL,    &p_nopbc,         0, "Ignore pbcs",                 0},
+    {"blk_avg",      '\0', POPT_ARG_NONE | POPT_ARGFLAG_OPTIONAL,   &p_blk_avg,       0, "Block Average",               0},
+    POPT_TABLEEND
+  };
+
+  // parse options
+  poptContext opt_context = poptGetContext("resdata", argc, argv, optionsTable, 0);
+  int opt=poptGetNextOpt(opt_context); // needs to be run to parse
+  if (opt < -1) {
+      /* Handle error condition */
+      fprintf(stderr, "%s: %s\n", poptBadOption(opt_context, POPT_BADOPTION_NOALIAS), poptStrerror(opt));
+      return 1;
+  }
+  poptFreeContext(opt_context);
+
+  // check if traj and top are set
+  if ( !(p_traj_path && p_top_path) )
+  {
+    std::cerr << "Trajectory and topology files must be set!" << std::endl;
+    return 1;
+  }
+
+  traj_path = std::string(p_traj_path);
+  top_path = std::string(p_top_path);
+  if ( p_index_path != NULL ) index_path = std::string(p_index_path);
+  mode = p_mode ? std::string(p_mode) : std::string("intra+same+cross");
+  if ( p_weights_path != NULL ) weights_path = std::string(p_weights_path);
+  if ( p_out_prefix != NULL ) out_prefix = std::string(p_out_prefix);
+  if ( p_nopbc != NULL ) nopbc = true;
+  if ( p_blk_avg != NULL ) blk_avg = true;
+  // check if paths are valid
+  if ( !std::filesystem::exists(std::filesystem::path(traj_path)) )
+  {
+    std::cerr << "Trajectory file does not exist!" << std::endl;
+    return 1;
+  }
+  if ( !std::filesystem::exists(std::filesystem::path(top_path)) )
+  {
+    std::cerr << "Topology file does not exist!" << std::endl;
+    return 2;
+  }
+  if ( !weights_path.empty() && !std::filesystem::exists(std::filesystem::path(weights_path)) )
+  {
+    std::cerr << "Weights file does not exist!" << std::endl;
+    return 3;
+  }
+  if (  !index_path.empty() &&!std::filesystem::exists(std::filesystem::path(index_path)) )
+  {
+    std::cerr << "Index file does not exist!" << std::endl;
+    return 4;
+  }
+  // if ( !out_prefix.empty() )
+  // {
+  //   bool created = std::filesystem::create_directories(std::filesystem::path(out_prefix));
+  //   if ( !created ) // if not created
+  //   {
+  //     std::cout << "Could not create output directory at " << out_prefix << std::endl;
+  //     if ( std::filesystem::exists(std::filesystem::path(out_prefix)) ) // already exists
+  //     {
+  //       std::cout << "Reason: directory already exists! WARNING: Files might be overwritten!" << std::endl;
+  //     }
+  //     else if ( !std::filesystem::is_directory(std::filesystem::path(out_prefix)) ) // not a directory (file or non-existent path)
+  //     {
+  //       std::cout << "Reason: path is not a directory!" << std::endl;
+  //     }
+  //     else if ( !std::filesystem::exists(std::filesystem::path(out_prefix)) ) // does not exist (no permissions or non-existent parent directory)
+  //     {
+  //       std::cout << "Reason: could not create directory!" << std::endl;
+  //       return 5;
+  //     }
+  //   }
+  // }
+  if ( num_threads < 1 )
+  {
+    std::cerr << "Number of threads must be at least 1!" << std::endl;
+    return 6;
+  }
+  if ( mol_threads < 1 )
+  {
+    std::cout << "Setting molecule threads to number of threads!" << std::endl;
+    mol_threads = num_threads;
+  }
+  if ( dt < 0 )
+  {
+    std::cerr << "Time step must be a positive number!" << std::endl;
+    return 7;
+  }
+  if ( nskip < 0 )
+  {
+    std::cerr << "Number of frames to skip must be at least 0!" << std::endl;
+    return 8;
+  }
+  if ( cutoff <= 0.0 )
+  {
+    std::cerr << "Cutoff distance must be greater than 0!" << std::endl;
+    return 9;
+  }
+  if ( mol_cutoff <= 0.0 )
+  {
+    std::cerr << "Molecule cutoff distance must be greater than 0!" << std::endl;
+    return 10;
+  }
+  if ( t_begin < 0.0 )
+  {
+    std::cerr << "Start time must be at least 0!" << std::endl;
+    return 11;
+  }
+  if ( t_end < t_begin && t_end != -1.0 )
+  {
+    std::cerr << "End time must be greater than start time!" << std::endl;
+    return 12;
+  }
+
+  resdata::RESData resdata(
+    top_path, traj_path, cutoff, mol_cutoff, nskip, num_threads, mol_threads, dt,
+    mode, weights_path, nopbc, t_begin, t_end, index_path, blk_avg
+  );
+  resdata.run();
+  // cmdata.process_data();
+  resdata.write_output(out_prefix);
+
+  return 0;
+}

tools/resdata/popt.patch

@@ -0,0 +1,42 @@
+diff --color -Naur CMakeLists.txt.origCMakeLists.txt
+--- CMakeLists.txt	2024-03-14 16:24:38.310420219 +0100
++++ CMakeLists.txt	2024-03-14 16:25:42.162353238 +0100
+@@ -103,12 +103,6 @@
+
+
+ add_subdirectory(src)
+-add_subdirectory(doc)
+ if (EXISTS po/popt.pot)
+ 	add_subdirectory(po)
+ endif()
+-
+-# Enable testing
+-include(CTest)
+-enable_testing()
+-add_subdirectory(tests)
+diff --color -Naur src/CMakeLists.txt.orig src/CMakeLists.txt
+--- src/CMakeLists.txt	2024-03-14 16:23:52.010463824 +0100
++++ src/CMakeLists.txt	2024-03-14 16:26:29.646298783 +0100
+@@ -25,13 +25,15 @@
+ endif()
+
+ set_target_properties(popt PROPERTIES
+-	VERSION ${PROJECT_VERSION}
+-	SOVERSION ${POPT_SOVERSION}
+-	C_STANDARD 99
+-	C_STANDARD_REQUIRED ON
+-	C_EXTENSIONS ON
+-	PUBLIC_HEADER popt.h
+-	LINK_FLAGS "-Wl,--no-undefined -Wl,--version-script,\"${PROJECT_SOURCE_DIR}/src/libpopt.vers\""
++    VERSION ${PROJECT_VERSION}
++    SOVERSION ${POPT_SOVERSION}
++    C_STANDARD 99
++    C_STANDARD_REQUIRED ON
++    C_EXTENSIONS ON
++    PUBLIC_HEADER popt.h
++    if(not APPLE)
++        LINK_FLAGS "-Wl,--no-undefined -Wl,--version-script,\"${PROJECT_SOURCE_DIR}/src/libpopt.vers\""
++    endif()
+ )
+
+ # Install the library