diff --git a/.github/docker/Dockerfile.centos-stream9 b/.github/docker/Dockerfile.centos-stream9 index f60341c..2389a28 100644 --- a/.github/docker/Dockerfile.centos-stream9 +++ b/.github/docker/Dockerfile.centos-stream9 @@ -1,4 +1,4 @@ FROM quay.io/centos/centos:stream9 AS base # Install dependencies -RUN yum install -y popt-devel gcc-c++ make diffutils +RUN yum install -y gcc-c++ make diffutils diff --git a/.github/docker/Dockerfile.debian-10 b/.github/docker/Dockerfile.debian-10 index 864ca86..e754f51 100644 --- a/.github/docker/Dockerfile.debian-10 +++ b/.github/docker/Dockerfile.debian-10 @@ -1,4 +1,4 @@ FROM debian:10 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ diff --git a/.github/docker/Dockerfile.debian-11 b/.github/docker/Dockerfile.debian-11 index 35cba82..8501ac9 100644 --- a/.github/docker/Dockerfile.debian-11 +++ b/.github/docker/Dockerfile.debian-11 @@ -1,4 +1,4 @@ FROM debian:11 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ diff --git a/.github/docker/Dockerfile.debian-12 b/.github/docker/Dockerfile.debian-12 index b8b09e5..678e303 100644 --- a/.github/docker/Dockerfile.debian-12 +++ b/.github/docker/Dockerfile.debian-12 @@ -1,4 +1,4 @@ FROM debian:12 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ diff --git a/.github/docker/Dockerfile.ubuntu-20.04 b/.github/docker/Dockerfile.ubuntu-20.04 index 3785af8..7fccca5 100644 --- a/.github/docker/Dockerfile.ubuntu-20.04 +++ b/.github/docker/Dockerfile.ubuntu-20.04 @@ -1,4 +1,4 @@ FROM ubuntu:20.04 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ diff --git a/.github/docker/Dockerfile.ubuntu-22.04 b/.github/docker/Dockerfile.ubuntu-22.04 index eaa4c2a..6f2b116 100644 --- a/.github/docker/Dockerfile.ubuntu-22.04 +++ b/.github/docker/Dockerfile.ubuntu-22.04 @@ -1,4 +1,4 @@ FROM ubuntu:24.04 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ \ No newline at end of file diff --git a/.github/docker/Dockerfile.ubuntu-22.04-clang b/.github/docker/Dockerfile.ubuntu-22.04-clang index 78c0616..1e4289b 100644 --- a/.github/docker/Dockerfile.ubuntu-22.04-clang +++ b/.github/docker/Dockerfile.ubuntu-22.04-clang @@ -1,4 +1,4 @@ FROM ubuntu:22.04 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make clang libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make clang diff --git a/.github/docker/Dockerfile.ubuntu-24.04 b/.github/docker/Dockerfile.ubuntu-24.04 index bd411c9..26c215b 100644 --- a/.github/docker/Dockerfile.ubuntu-24.04 +++ b/.github/docker/Dockerfile.ubuntu-24.04 @@ -1,4 +1,4 @@ FROM ubuntu:22.04 AS base # Install dependencies -RUN apt-get update && apt-get install -y --no-install-recommends make g++ libpopt-dev libpopt0 +RUN apt-get update && apt-get install -y --no-install-recommends make g++ diff --git a/Makefile b/Makefile index 1f48364..43e25e6 100644 --- a/Makefile +++ b/Makefile @@ -95,11 +95,11 @@ endif DEFINES += -DRBT_VERSION=\"$(RBT_VERSION)\" CXX_FLAGS := $(CXX_BASE_FLAGS) $(CXX_CONFIG_FLAGS) $(CXX_WARNING_FLAGS) $(CXX_EXTRA_FLAGS) $(DEFINES) LINK_FLAGS := -shared -LIB_DEPENDENCIES = -lm -lpopt +LIB_DEPENDENCIES += -lm # Modify library dependencies for macOS ifeq ($(shell uname),Darwin) - LIB_DEPENDENCIES = -L$(BREW_PREFIX)/lib -lpopt -lm + LIB_DEPENDENCIES = -L$(BREW_PREFIX)/lib -lm endif LIBS += $(LIB_DEPENDENCIES) -lRbt diff --git a/README.md b/README.md index 43e5948..55a0950 100644 --- a/README.md +++ b/README.md @@ -22,13 +22,12 @@ make sure the following requirements are installed and available: - make - a c++ compiler (g++ by default) -- popt and development headers (libpopt0 and libpopt-dev in ubuntu) - git (optional if you download the code directly) if you're running ubuntu, you can get all of them by running ``` -sudo apt update && sudo apt install -y make git libpopt0 libpopt-dev g++ +sudo apt update && sudo apt install -y make git g++ ``` if you're running `macOS`, make sure to use gcc instead of clang. You can install gcc and popt using homebrew diff --git a/include/RbtArgParser.h b/include/RbtArgParser.h index 285b118..d85912b 100644 --- a/include/RbtArgParser.h +++ b/include/RbtArgParser.h @@ -2,6 +2,7 @@ #define _RBTARGPARSER_H_ #include +#include #include #include @@ -61,6 +62,12 @@ class RbtOptionValue { return *this; } + template + RbtOptionValue &operator>>(std::optional &v) { + v = value.as(); + return *this; + } + inline bool is_present() { return value.count() > 0; } }; diff --git a/include/rbdock/rbdock_argparser.h b/include/rbdock/rbdock_argparser.h new file mode 100644 index 0000000..5dd23f4 --- /dev/null +++ b/include/rbdock/rbdock_argparser.h @@ -0,0 +1,16 @@ +#ifndef _RBT_RBDOCK_ARGPARSER_H_ +#define _RBT_RBDOCK_ARGPARSER_H_ + +#include "rbdock/rbdock_config.h" + +#include "RbtArgParser.h" + +namespace RBDock { + +RbtArgParser::RbtArgParser get_options_parser(); + +RBDock::RBDockConfig parse_args(int argc, const char *argv[]); + +} // namespace RBDock + +#endif // _RBT_RBDOCK_ARGPARSER_H_ \ No newline at end of file diff --git a/include/rbdock/rbdock_config.h b/include/rbdock/rbdock_config.h new file mode 100644 index 0000000..916827e --- /dev/null +++ b/include/rbdock/rbdock_config.h @@ -0,0 +1,114 @@ +#ifndef _RBDOCK_CONFIG_H_ +#define _RBDOCK_CONFIG_H_ + +#include +#include +#include +#include + +#include "RbtValidationError.h" + +namespace RBDock { + +struct RBDockConfig { + // mandatory parameters + std::string strLigandMdlFile; + std::string strReceptorPrmFile; + std::string strParamFile; + + // optional parameters + std::optional strOutputPrefix; + std::optional strFilterFile; + std::optional nDockingRuns; + std::optional dTargetScore; + std::optional nSeed; + std::optional iTrace; + + // flags + bool bContinue = false; + bool bPosIonise = false; + bool bNegIonise = false; + bool bAllH = false; + + friend std::ostream &operator<<(std::ostream &os, const RBDockConfig &config) { + os << "Command line args:" << std::endl; + os << " -i " << config.strLigandMdlFile << std::endl; + os << " -r " << config.strReceptorPrmFile << std::endl; + os << " -p " << config.strParamFile << std::endl; + if (config.strOutputPrefix.has_value()) { + os << " -o " << config.strOutputPrefix.value() << std::endl; + } else { + os << "WARNING: output file name is missing." << std::endl; + } + auto def_aux = config.nDockingRuns.has_value() ? " (default) " : ""; + os << " -n " << config.nDockingRuns.value_or(1) << def_aux << std::endl; + if (config.nSeed.has_value()) os << " -s " << config.nSeed.value() << std::endl; + if (config.iTrace.has_value()) os << " -T " << config.iTrace.value() << std::endl; + if (config.bPosIonise) os << " --ap " << std::endl; + if (config.bNegIonise) os << " --an " << std::endl; + if (!config.bAllH) os << " --allH " << std::endl; + if (config.dTargetScore.has_value()) os << " -t " << config.dTargetScore.value() << std::endl; + if (config.bContinue) os << " --cont " << std::endl; + if (config.bPosIonise) + os << "Automatically protonating positive ionisable groups (amines, imidazoles, guanidines)" << std::endl; + if (config.bNegIonise) + os << "Automatically deprotonating negative ionisable groups (carboxylic acids, phosphates, sulphates, " + "sulphonates)" + << std::endl; + if (!config.bAllH) + os << "Reading polar hydrogens only from ligand SD file" << std::endl; + else + os << "Reading all hydrogens from ligand SD file" << std::endl; + + if (config.dTargetScore.has_value()) { + os << std::endl << "Lower target intermolecular score = " << config.dTargetScore.value() << std::endl; + } + if (config.strFilterFile.has_value()) os << " -f " << config.strFilterFile.value() << std::endl; + + return os; + } + + void validate() { + if (strLigandMdlFile.empty()) throw Rbt::ValidationError("input ligand SD file is mandatory"); + if (strReceptorPrmFile.empty()) throw Rbt::ValidationError("receptor parameter file is mandatory"); + if (nDockingRuns.has_value() && nDockingRuns < 1) + throw Rbt::ValidationError("number of runs must be greater than 0"); + } + + bool operator==(const RBDockConfig &rhs) const { + return strLigandMdlFile == rhs.strLigandMdlFile && strOutputPrefix == rhs.strOutputPrefix + && strReceptorPrmFile == rhs.strReceptorPrmFile && strParamFile == rhs.strParamFile + && strFilterFile == rhs.strFilterFile && nDockingRuns == rhs.nDockingRuns && bContinue == rhs.bContinue + && dTargetScore == rhs.dTargetScore && bPosIonise == rhs.bPosIonise && bNegIonise == rhs.bNegIonise + && bAllH == rhs.bAllH && nSeed == rhs.nSeed && iTrace == rhs.iTrace; + } + + std::string get_filter_string() const { + std::ostringstream strAuxFilter; + if (strFilterFile.has_value()) return ""; + if (dTargetScore.has_value()) { // -t + if (!nDockingRuns.has_value()) { // -t only + strAuxFilter << "0 1 - SCORE.INTER " << dTargetScore.value() << std::endl; + } else // -t -n need to check if -cont present + // for all other cases it doesn't matter + if (bContinue) { // -t -n -cont + strAuxFilter << "1 if - SCORE.NRUNS " << (nDockingRuns.value() - 1) + << " 0.0 -1.0,\n1 - SCORE.INTER " << dTargetScore.value() << std::endl; + } else { // -t -n + strAuxFilter << "1 if - " << dTargetScore.value() << " SCORE.INTER 0.0 " + << "if - SCORE.NRUNS " << (nDockingRuns.value() - 1) << " 0.0 -1.0,\n1 - SCORE.INTER " + << dTargetScore.value() << std::endl; + } + } // no target score, no filter + else if (nDockingRuns.has_value()) { // -n + strAuxFilter << "1 if - SCORE.NRUNS " << (nDockingRuns.value() - 1) << " 0.0 -1.0,\n0"; + } else { // no -t no -n + strAuxFilter << "0 0\n"; + } + return strAuxFilter.str(); + } +}; + +} // namespace RBDock + +#endif // _RBDOCK_CONFIG_H_ diff --git a/include/rbdock/rbdock_main.h b/include/rbdock/rbdock_main.h new file mode 100644 index 0000000..49a3caf --- /dev/null +++ b/include/rbdock/rbdock_main.h @@ -0,0 +1,12 @@ +#ifndef _RBT_RBDOCK_MAIN_H_ +#define _RBT_RBDOCK_MAIN_H_ + +#include + +#include "rbdock/rbdock_config.h" + +namespace RBDock { +void RBDock(const RBDockConfig &config, const std::string &strExeName); +} // namespace RBDock + +#endif // _RBT_RBDOCK_MAIN_H_ \ No newline at end of file diff --git a/src/exe/rbdock.cxx b/src/exe/rbdock.cxx index e6ae323..698d5db 100644 --- a/src/exe/rbdock.cxx +++ b/src/exe/rbdock.cxx @@ -9,517 +9,41 @@ * the University of Barcelona. * http://rdock.sourceforge.net/ ***********************************************************************/ +#include "rbdock/rbdock_argparser.h" +#include "rbdock/rbdock_config.h" +#include "rbdock/rbdock_main.h" -// Main docking application -#include -using std::setw; - -#include -#include // for command-line parsing - -#include "RbtBiMolWorkSpace.h" -#include "RbtCrdFileSink.h" -#include "RbtDockingError.h" -#include "RbtFileError.h" -#include "RbtFilter.h" -#include "RbtLigandError.h" -#include "RbtMdlFileSink.h" -#include "RbtMdlFileSource.h" -#include "RbtModelError.h" -#include "RbtPRMFactory.h" -#include "RbtParameterFileSource.h" -#include "RbtRand.h" -#include "RbtSFFactory.h" -#include "RbtSFRequest.h" -#include "RbtTransformFactory.h" +#include "Rbt.h" +#include "RbtDebug.h" #include "RbtVersion.h" const RbtString EXEVERSION = RBT_VERSION; -// Section name in docking prm file containing scoring function definition -const RbtString _ROOT_SF = "SCORE"; -const RbtString _RESTRAINT_SF = "RESTR"; -const RbtString _ROOT_TRANSFORM = "DOCK"; - -void PrintUsage(void) { - cout << endl << "Usage:" << endl; - cout << "rbdock -i -o -r -p [-n ] [-ap] [-an] [-allH]" - << endl; - cout << " [-t ] [-c] [-T ] [-s ]" << endl; - cout << endl << "Options:\t-i - input ligand SD file" << endl; - cout << "\t\t-o - root name for output file(s)" << endl; - cout << "\t\t-r - receptor parameter file " << endl; - cout << "\t\t-p - docking protocol parameter file" << endl; - cout << "\t\t-n - number of runs/ligand (default=1)" << endl; - cout << "\t\t-ap - protonate all neutral amines, guanidines, imidazoles (default=disabled)" << endl; - cout << "\t\t-an - deprotonate all carboxylic, sulphur and phosphorous acid groups (default=disabled)" << endl; - cout << "\t\t-allH - read all hydrogens present (default=polar hydrogens only)" << endl; - cout << "\t\t-t - score threshold OR filter file name" << endl; - cout << "\t\t-c - continue if score threshold is met (use with -t , default=terminate ligand)" - << endl; - cout << "\t\t-T - controls output level for debugging (0 = minimal, >0 = more verbose)" << endl; - cout << "\t\t-s - random number seed (default=from sys clock)" << endl; -} ///////////////////////////////////////////////////////////////////// // MAIN PROGRAM STARTS HERE ///////////////////////////////////////////////////////////////////// int main(int argc, const char *argv[]) { - cout.setf(ios_base::left, ios_base::adjustfield); - // Strip off the path to the executable, leaving just the file name - RbtString strExeName(argv[0]); - RbtString::size_type i = strExeName.rfind("/"); - if (i != RbtString::npos) strExeName.erase(0, i + 1); - - // Print a standard header - Rbt::PrintStdHeader(cout, strExeName + " - " + EXEVERSION); - - // Command line arguments and default values - RbtString strLigandMdlFile; - RbtBool bOutput(false); - RbtString strRunName; - RbtString strReceptorPrmFile; // Receptor param file - RbtString strParamFile; // Docking run param file - RbtString strFilterFile; // Filter file - RbtInt nDockingRuns(0); // Init to zero, so can detect later whether user explictly typed -n - - // Params for target score - RbtBool bTarget(false); - RbtBool bStop(true); // DM 25 May 2001 - if true, stop once target is met - RbtBool bDockingRuns(false); // is argument -n present? - RbtDouble dTargetScore(0.0); - RbtBool bFilter(false); - - RbtBool bPosIonise(false); - RbtBool bNegIonise(false); - RbtBool bImplH(true); // if true, read only polar hydrogens from SD file, else read all H's present - RbtBool bSeed(false); // Random number seed (default = from system clock) - RbtInt nSeed(0); - RbtBool bTrace(false); - RbtInt iTrace(0); // Trace level, for debugging - - // variables for popt command-line parsing - char c; // for argument parsing - poptContext optCon; // ditto - - char *inputFile = NULL; // will be 'strLigandMdlFile' - char *outputFile = NULL; // will be 'strRunName' - char *receptorFile = NULL; // will be 'strReceptorPrmFile' - char *protocolFile = NULL; // will be 'strParamFile' - char *strTargetScr = NULL; // will be 'dTargetScore' - struct poptOption optionsTable[] = { - // command line options - {"input", 'i', POPT_ARG_STRING | POPT_ARGFLAG_ONEDASH, &inputFile, 'i', "input file"}, - {"output", 'o', POPT_ARG_STRING | POPT_ARGFLAG_ONEDASH, &outputFile, 'o', "output file"}, - {"receptor", 'r', POPT_ARG_STRING | POPT_ARGFLAG_ONEDASH, &receptorFile, 'r', "receptor file"}, - {"protocol", 'p', POPT_ARG_STRING | POPT_ARGFLAG_ONEDASH, &protocolFile, 'p', "protocol file"}, - {"runs", 'n', POPT_ARG_INT | POPT_ARGFLAG_ONEDASH, &nDockingRuns, 'n', "number of runs"}, - {"trace", 'T', POPT_ARG_INT | POPT_ARGFLAG_ONEDASH, &iTrace, 'T', "trace level for debugging"}, - {"seed", 's', POPT_ARG_INT | POPT_ARGFLAG_ONEDASH, &nSeed, 's', "random seed"}, - {"ap", 'P', POPT_ARG_NONE | POPT_ARGFLAG_ONEDASH, 0, 'P', "protonate groups"}, - {"an", 'D', POPT_ARG_NONE | POPT_ARGFLAG_ONEDASH, 0, 'D', "DEprotonate groups"}, - {"allH", 'H', POPT_ARG_NONE | POPT_ARGFLAG_ONEDASH, 0, 'H', "read all Hs"}, - {"target", 't', POPT_ARG_STRING | POPT_ARGFLAG_ONEDASH, &strTargetScr, 't', "target score"}, - {"cont", 'C', POPT_ARG_NONE | POPT_ARGFLAG_ONEDASH, 0, 'C', "continue even if target met"}, - POPT_AUTOHELP{NULL, 0, 0, NULL, 0}}; - - optCon = poptGetContext(NULL, argc, argv, optionsTable, 0); - poptSetOtherOptionHelp(optCon, "-r -p -i -o [options]"); - // Brief help message - if (argc < 2) { - PrintUsage(); - return 1; - } - RbtDouble val(0.0); - while ((c = poptGetNextOpt(optCon)) >= 0) { - switch (c) { - case 'P': // protonate - bPosIonise = true; - break; - case 'D': // deprotonate - bNegIonise = true; - break; - case 'H': // all-H model - bImplH = false; - break; - case 'C': - bStop = false; - break; - case 't': - // If str can be translated to an integer, I assume is a - // threshold. Otherwise, I assume is the filter file name - char *error; - errno = 0; - val = strtod(strTargetScr, &error); - if (!errno && !*error) // Is it a number? - { - dTargetScore = val; - bTarget = true; - } else // Assume is the filter file name - { - bFilter = true; - strFilterFile = strTargetScr; - } - break; - case 's': - bSeed = true; - break; - case 'T': - bTrace = true; - break; - default: - break; - } - } - cout << endl; - poptFreeContext(optCon); - - // print out arguments - // input ligand file, receptor and parameter is compulsory - cout << endl << "Command line args:" << endl; - if (!inputFile || !receptorFile || !protocolFile) { // if any of them is missing - // SRC poptPrintUsage(optCon, stderr, 0); - PrintUsage(); - exit(1); - } else { - strLigandMdlFile = inputFile; - strReceptorPrmFile = receptorFile; - strParamFile = protocolFile; - cout << " -i " << strLigandMdlFile << endl; - cout << " -r " << strReceptorPrmFile << endl; - cout << " -p " << strParamFile << endl; - } - // output is not that important but good to have - if (outputFile) { - strRunName = outputFile; - bOutput = true; - cout << " -o " << strRunName << endl; - } else { - cout << "WARNING: output file name is missing." << endl; - } - // docking runs - if (nDockingRuns >= 1) { // User typed -n explicitly, so set bDockingRuns to true - bDockingRuns = true; - cout << " -n " << nDockingRuns << endl; - } else { - nDockingRuns = 1; // User didn't type -n explicitly, so fall back to the default of n=1 - cout << " -n " << nDockingRuns << " (default) " << endl; - } - if (bSeed) // random seed (if provided) - cout << " -s " << nSeed << endl; - if (bTrace) // random seed (if provided) - cout << " -T " << iTrace << endl; - if (bPosIonise) // protonate - cout << " -ap " << endl; - if (bNegIonise) // deprotonate - cout << " -an " << endl; - if (!bImplH) // all-H - cout << " -allH " << endl; - if (!bStop) // stop after target - cout << " -cont " << endl; - if (bTarget) cout << " -t " << dTargetScore << endl; - - // BGD 26 Feb 2003 - Create filters to simulate old rbdock - // behaviour - ostringstream strFilter; - if (!bFilter) { - if (bTarget) // -t - { - if (!bDockingRuns) // -t only - { - strFilter << "0 1 - SCORE.INTER " << dTargetScore << endl; - } else // -t -n need to check if -cont present - // for all other cases it doesn't matter - if (!bStop) // -t -n -cont - { - strFilter << "1 if - SCORE.NRUNS " << (nDockingRuns - 1) << " 0.0 -1.0,\n1 - SCORE.INTER " - << dTargetScore << endl; - } else // -t -n - { - strFilter << "1 if - " << dTargetScore << " SCORE.INTER 0.0 " - << "if - SCORE.NRUNS " << (nDockingRuns - 1) << " 0.0 -1.0,\n1 - SCORE.INTER " - << dTargetScore << endl; - } - } // no target score, no filter - else if (bDockingRuns) // -n - { - strFilter << "1 if - SCORE.NRUNS " << (nDockingRuns - 1) << " 0.0 -1.0,\n0"; - } else // no -t no -n - { - strFilter << "0 0\n"; - } - } - - // DM 20 Apr 1999 - set the auto-ionise flags - if (bPosIonise) - cout << "Automatically protonating positive ionisable groups (amines, imidazoles, guanidines)" << endl; - if (bNegIonise) - cout << "Automatically deprotonating negative ionisable groups (carboxylic acids, phosphates, sulphates, " - "sulphonates)" - << endl; - if (bImplH) - cout << "Reading polar hydrogens only from ligand SD file" << endl; - else - cout << "Reading all hydrogens from ligand SD file" << endl; - - if (bTarget) { - cout << endl << "Lower target intermolecular score = " << dTargetScore << endl; - } + std::string exeFullPath(argv[0]); + std::string strExeName = exeFullPath.substr(exeFullPath.find_last_of("/\\") + 1); + Rbt::PrintStdHeader(std::cout, strExeName + " - " + EXEVERSION); try { - // Create a bimolecular workspace - RbtBiMolWorkSpacePtr spWS(new RbtBiMolWorkSpace()); - // Set the workspace name to the root of the receptor .prm filename - RbtStringList componentList = Rbt::ConvertDelimitedStringToList(strReceptorPrmFile, "."); - RbtString wsName = componentList.front(); - spWS->SetName(wsName); - - // Read the docking protocol parameter file - RbtParameterFileSourcePtr spParamSource( - new RbtParameterFileSource(Rbt::GetRbtFileName("data/scripts", strParamFile)) - ); - // Read the receptor parameter file - RbtParameterFileSourcePtr spRecepPrmSource( - new RbtParameterFileSource(Rbt::GetRbtFileName("data/receptors", strReceptorPrmFile)) - ); - cout << endl - << "DOCKING PROTOCOL:" << endl - << spParamSource->GetFileName() << endl - << spParamSource->GetTitle() << endl; - cout << endl - << "RECEPTOR:" << endl - << spRecepPrmSource->GetFileName() << endl - << spRecepPrmSource->GetTitle() << endl; - - // Create the scoring function from the SCORE section of the docking protocol prm file - // Format is: - // SECTION SCORE - // INTER RbtInterSF.prm - // INTRA RbtIntraSF.prm - // END_SECTION - // - // Notes: - // Section name must be SCORE. This is also the name of the root SF aggregate - // An aggregate is created for each parameter in the section. - // Parameter name becomes the name of the subaggregate (e.g. SCORE.INTER) - // Parameter value is the file name for the subaggregate definition - // Default directory is $RBT_ROOT/data/sf - RbtSFFactoryPtr spSFFactory(new RbtSFFactory()); // Factory class for scoring functions - RbtSFAggPtr spSF(new RbtSFAgg(_ROOT_SF)); // Root SF aggregate - spParamSource->SetSection(_ROOT_SF); - RbtStringList sfList(spParamSource->GetParameterList()); - // Loop over all parameters in the SCORE section - for (RbtStringListConstIter sfIter = sfList.begin(); sfIter != sfList.end(); sfIter++) { - // sfFile = file name for scoring function subaggregate - RbtString sfFile(Rbt::GetRbtFileName("data/sf", spParamSource->GetParameterValueAsString(*sfIter))); - RbtParameterFileSourcePtr spSFSource(new RbtParameterFileSource(sfFile)); - // Create and add the subaggregate - spSF->Add(spSFFactory->CreateAggFromFile(spSFSource, *sfIter)); - } - - // Add the RESTRAINT subaggregate scoring function from any SF definitions in the receptor prm file - spSF->Add(spSFFactory->CreateAggFromFile(spRecepPrmSource, _RESTRAINT_SF)); - - // Create the docking transform aggregate from the transform definitions in the docking prm file - RbtTransformFactoryPtr spTransformFactory(new RbtTransformFactory()); - spParamSource->SetSection(); - RbtTransformAggPtr spTransform(spTransformFactory->CreateAggFromFile(spParamSource, _ROOT_TRANSFORM)); - - // Override the TRACE levels for the scoring function and transform - // Dump details to cout - // Register the scoring function and the transform with the workspace - if (bTrace) { - RbtRequestPtr spTraceReq(new RbtSFSetParamRequest("TRACE", iTrace)); - spSF->HandleRequest(spTraceReq); - spTransform->HandleRequest(spTraceReq); - } - if (iTrace > 0) { - cout << endl << "SCORING FUNCTION DETAILS:" << endl << *spSF << endl; - cout << endl << "SEARCH DETAILS:" << endl << *spTransform << endl; - } - spWS->SetSF(spSF); - spWS->SetTransform(spTransform); - - // DM 18 May 1999 - // Variants describing the library version, exe version, parameter file, and current directory - // Will be stored in the ligand SD files - RbtVariant vLib(Rbt::GetProduct() + " (" + Rbt::GetVersion() + ", Build" + Rbt::GetBuild() + ")"); - RbtVariant vExe(strExeName + " - " + EXEVERSION); - RbtVariant vRecep(spRecepPrmSource->GetFileName()); - RbtVariant vPrm(spParamSource->GetFileName()); - RbtVariant vDir(Rbt::GetCurrentDirectory()); - - spRecepPrmSource->SetSection(); - // Read docking site from file and register with workspace - RbtString strASFile = spWS->GetName() + ".as"; - RbtString strInputFile = Rbt::GetRbtFileName("data/grids", strASFile); - // DM 26 Sep 2000 - ios_base::binary is invalid with IRIX CC -#if defined(__sgi) && !defined(__GNUC__) - ifstream istr(strInputFile.c_str(), ios_base::in); -#else - ifstream istr(strInputFile.c_str(), ios_base::in | ios_base::binary); -#endif - // DM 14 June 2006 - bug fix to one of the longest standing rDock issues - //(the cryptic "Error reading from input stream" message, if cavity file was missing) - if (!istr) { - RbtString message = "Cavity file (" + strASFile + ") not found in current directory or $RBT_HOME"; - message += " - run rbcavity first"; - throw RbtFileReadError(_WHERE_, message); - } - RbtDockingSitePtr spDS(new RbtDockingSite(istr)); - istr.close(); - spWS->SetDockingSite(spDS); - cout << endl << "DOCKING SITE" << endl << (*spDS) << endl; - - // Prepare the SD file sink for saving the docked conformations for each ligand - // DM 3 Dec 1999 - replaced ostringstream with RbtString in determining SD file name - // SRC 2014 moved here this block to allow WRITE_ERROR TRUE - if (bOutput) { - RbtMolecularFileSinkPtr spMdlFileSink(new RbtMdlFileSink(strRunName + ".sd", RbtModelPtr())); - spWS->SetSink(spMdlFileSink); - } - - RbtPRMFactory prmFactory(spRecepPrmSource, spDS); - prmFactory.SetTrace(iTrace); - // Create the receptor model from the file names in the receptor parameter file - RbtModelPtr spReceptor = prmFactory.CreateReceptor(); - spWS->SetReceptor(spReceptor); - - // Register any solvent - RbtModelList solventList = prmFactory.CreateSolvent(); - spWS->SetSolvent(solventList); - if (spWS->hasSolvent()) { - RbtInt nSolvent = spWS->GetSolvent().size(); - cout << endl << nSolvent << " solvent molecules registered" << endl; - } else { - cout << endl << "No solvent" << endl; - } - - // SRC 2014 removed sector bOutput from here to some blocks above, for WRITEERRORS TRUE - - // Seed the random number generator - RbtRand &theRand = Rbt::GetRbtRand(); // ref to random number generator - if (bSeed) { - theRand.Seed(nSeed); - } - - // Create the filter object for controlling early termination of protocol - RbtFilterPtr spfilter; - if (bFilter) { - spfilter = new RbtFilter(strFilterFile); - if (bDockingRuns) { - spfilter->SetMaxNRuns(nDockingRuns); - } - } else { - spfilter = new RbtFilter(strFilter.str(), true); - } - if (bTrace) { - RbtRequestPtr spTraceReq(new RbtSFSetParamRequest("TRACE", iTrace)); - spfilter->HandleRequest(spTraceReq); - } - - // Register the Filter with the workspace - spWS->SetFilter(spfilter); - - // MAIN LOOP OVER LIGAND RECORDS - // DM 20 Apr 1999 - add explicit bPosIonise and bNegIonise flags to MdlFileSource constructor - RbtMolecularFileSourcePtr spMdlFileSource( - new RbtMdlFileSource(strLigandMdlFile, bPosIonise, bNegIonise, bImplH) - ); - for (RbtInt nRec = 1; spMdlFileSource->FileStatusOK(); spMdlFileSource->NextRecord(), nRec++) { - cout.setf(ios_base::left, ios_base::adjustfield); - cout << endl << "**************************************************" << endl << "RECORD #" << nRec << endl; - RbtError molStatus = spMdlFileSource->Status(); - if (!molStatus.isOK()) { - cout << endl << molStatus << endl << "************************************************" << endl; - continue; - } - - // DM 26 Jul 1999 - only read the largest segment (guaranteed to be called H) - // BGD 07 Oct 2002 - catching errors created by the ligands, - // so rbdock continues with the next one, instead of - // completely stopping - try { - spMdlFileSource->SetSegmentFilterMap(Rbt::ConvertStringToSegmentMap("H")); - - if (spMdlFileSource->isDataFieldPresent("Name")) - cout << "NAME: " << spMdlFileSource->GetDataValue("Name") << endl; - if (spMdlFileSource->isDataFieldPresent("REG_Number")) - cout << "REG_Num:" << spMdlFileSource->GetDataValue("REG_Number") << endl; - cout << setw(30) << "RANDOM_NUMBER_SEED:" << theRand.GetSeed() << endl; - - // Create and register the ligand model - RbtModelPtr spLigand = prmFactory.CreateLigand(spMdlFileSource); - RbtString strMolName = spLigand->GetName(); - spWS->SetLigand(spLigand); - // Update any model coords from embedded chromosomes in the ligand file - spWS->UpdateModelCoordsFromChromRecords(spMdlFileSource, iTrace); - - // DM 18 May 1999 - store run info in model data - // Clear any previous Rbt.* data fields - spLigand->ClearAllDataFields("Rbt."); - spLigand->SetDataValue("Rbt.Library", vLib); - spLigand->SetDataValue("Rbt.Executable", vExe); - spLigand->SetDataValue("Rbt.Receptor", vRecep); - spLigand->SetDataValue("Rbt.Parameter_File", vPrm); - spLigand->SetDataValue("Rbt.Current_Directory", vDir); - - // DM 10 Dec 1999 - if in target mode, loop until target score is reached - RbtBool bTargetMet = false; - - //////////////////////////////////////////////////// - // MAIN LOOP OVER EACH SIMULATED ANNEALING RUN - // Create a history file sink, just in case it's needed by any - // of the transforms - RbtInt iRun = 1; - // need to check this here. The termination - // filter is only run once at least - // one docking run has been done. - if (nDockingRuns < 1) bTargetMet = true; - while (!bTargetMet) { - // Catching errors with this specific run - try { - if (bOutput) { - ostringstream histr; - histr << strRunName << "_" << strMolName << nRec << "_his_" << iRun << ".sd"; - RbtMolecularFileSinkPtr spHistoryFileSink(new RbtMdlFileSink(histr.str(), spLigand)); - spWS->SetHistorySink(spHistoryFileSink); - } - spWS->Run(); // Dock! - RbtBool bterm = spfilter->Terminate(); - RbtBool bwrite = spfilter->Write(); - if (bterm) bTargetMet = true; - if (bOutput && bwrite) { - spWS->Save(); - } - iRun++; - } catch (RbtDockingError &e) { - cout << e << endl; - } - } - // END OF MAIN LOOP OVER EACH SIMULATED ANNEALING RUN - //////////////////////////////////////////////////// - } - // END OF TRY - catch (RbtLigandError &e) { - cout << e << endl; - } - } - // END OF MAIN LOOP OVER LIGAND RECORDS - //////////////////////////////////////////////////// - cout << endl << "END OF RUN" << endl; - // if (bOutput && flexRec) { - // RbtMolecularFileSinkPtr spRecepSink(new RbtCrdFileSink(strRunName+".crd",spReceptor)); - // spRecepSink->Render(); - // } + RBDock::RBDockConfig config = RBDock::parse_args(argc, argv); + std::cout << config << std::endl; + std::cout.setf(std::ios_base::left, std::ios_base::adjustfield); + RBDock::RBDock(config, strExeName); } catch (RbtError &e) { - cout << e << endl; - } catch (...) { - cout << "Unknown exception" << endl; + std::cerr << e << std::endl; + return 1; + } catch (std::exception &e) { + std::cerr << "Unknown exception" << std::endl; + std::cerr << typeid(e).name() << std::endl; + std::cerr << e.what() << std::endl; + return 1; } - _RBTOBJECTCOUNTER_DUMP_(cout) - + _RBTOBJECTCOUNTER_DUMP_(std::cout) return 0; } diff --git a/src/lib/rbdock/rbdock_argparser.cxx b/src/lib/rbdock/rbdock_argparser.cxx new file mode 100644 index 0000000..62b6a77 --- /dev/null +++ b/src/lib/rbdock/rbdock_argparser.cxx @@ -0,0 +1,76 @@ +#include "rbdock/rbdock_argparser.h" + +#include +#include + +#include "rbdock/rbdock_config.h" + +#include "RbtArgParser.h" +#include "RbtError.h" + +RbtArgParser::RbtArgParser RBDock::get_options_parser() { + using std::string; + + RbtArgParser::RbtArgParser parser("rbdock", "Docking application"); + parser.add("i,input", "input ligand SD file (mandatory)"); + parser.add("o,output-root", "root name for output file(s)", ""); + parser.add("r,receptor", "receptor parameter file (mandatory)"); + parser.add("p,protocol", "docking protocol parameter file", "dock.prm"); + parser.add("n,runs", "number of runs/ligand", "1"); + parser.add_flag("P,ap", "protonate all neutral amines, guanidines, imidazoles"); + parser.add_flag("D,an", "deprotonate all carboxylic, sulphur and phosphorous acid groups"); + parser.add_flag("H,allH", "read all hydrogens present. If disabled, read polar hydrogens only"); + parser.add("t,target", "score threshold (or filter file name, deprecated)"); + parser.add_flag("C,cont", "if enabled, continue if score threshold is met (use with -t )"); + parser.add("T,trace", "controls output level for debugging (0 = minimal, >0 = more verbose)", "0"); + parser.add("s,seed", "random number seed (default=from sys clock)"); + parser.add("f,filter", "filter file name", ""); + + return parser; +} + +RBDock::RBDockConfig RBDock::parse_args(int argc, const char *argv[]) { + auto parser = get_options_parser(); + try { + auto parser_result = RbtArgParser::RbtParseResult(parser.parse(argc, argv)); + RBDock::RBDockConfig config; + parser_result["input"] >> config.strLigandMdlFile; + parser_result["receptor"] >> config.strReceptorPrmFile; + if (parser_result["protocol"].is_present()) parser_result["protocol"] >> config.strParamFile; + if (parser_result["output-root"].is_present()) parser_result["output-root"] >> config.strOutputPrefix; + if (parser_result["runs"].is_present()) parser_result["runs"] >> config.nDockingRuns; + parser_result["ap"] >> config.bPosIonise; + parser_result["an"] >> config.bNegIonise; + parser_result["allH"] >> config.bAllH; + if (parser_result["seed"].is_present()) parser_result["seed"] >> config.nSeed; + if (parser_result["trace"].is_present()) parser_result["trace"] >> config.iTrace; + if (parser_result["target"].is_present()) { + std::string target; + parser_result["target"] >> target; + // if target is numeric, store as a target score + try { + config.dTargetScore = std::stod(target); + } catch (std::invalid_argument &e) { + // if target is not numeric, store as a filter file name + std::cerr << "Warning: -t option is deprecated for filter files. Use -f instead." << std::endl; + config.strFilterFile = target; + } + } + if (parser_result["filter"].is_present()) { + if (config.dTargetScore.has_value()) { + throw Rbt::ValidationError("Cannot specify both -t and -f options"); + } + parser_result["filter"] >> config.strFilterFile; + } + parser_result["cont"] >> config.bContinue; + + config.validate(); + return config; + } catch (RbtArgParser::ParsingError &e) { + std::cerr << "Error parsing options: " << e.what() << std::endl; + } catch (RbtArgParser::ValidationError &e) { + std::cerr << "Invalid configuration: " << e.what() << std::endl; + } + std::cerr << parser.help() << std::endl; + exit(1); +} \ No newline at end of file diff --git a/src/lib/rbdock/rbdock_main.cxx b/src/lib/rbdock/rbdock_main.cxx new file mode 100644 index 0000000..ac050c5 --- /dev/null +++ b/src/lib/rbdock/rbdock_main.cxx @@ -0,0 +1,277 @@ +#include "rbdock/rbdock_main.h" + +#include +#include + +#include "rbdock/rbdock_argparser.h" +#include "rbdock/rbdock_config.h" + +#include "RbtArgParser.h" +#include "RbtBiMolWorkSpace.h" +#include "RbtCrdFileSink.h" +#include "RbtDockingError.h" +#include "RbtFileError.h" +#include "RbtFilter.h" +#include "RbtLigandError.h" +#include "RbtMdlFileSink.h" +#include "RbtMdlFileSource.h" +#include "RbtModelError.h" +#include "RbtPRMFactory.h" +#include "RbtParameterFileSource.h" +#include "RbtPlatformCompatibility.h" +#include "RbtRand.h" +#include "RbtSFFactory.h" +#include "RbtSFRequest.h" +#include "RbtTransformFactory.h" +#include "RbtVersion.h" + +// Section name in docking prm file containing scoring function definition +const RbtString _ROOT_SF = "SCORE"; +const RbtString _RESTRAINT_SF = "RESTR"; +const RbtString _ROOT_TRANSFORM = "DOCK"; + +void RBDock::RBDock(const RBDock::RBDockConfig &config, const RbtString &strExeName) { + // Create a bimolecular workspace + RbtBiMolWorkSpacePtr spWS(new RbtBiMolWorkSpace()); + // Set the workspace name to the root of the receptor .prm filename + RbtStringList componentList = Rbt::ConvertDelimitedStringToList(config.strReceptorPrmFile, "."); + RbtString wsName = componentList.front(); + spWS->SetName(wsName); + + // Read the docking protocol parameter file + RbtParameterFileSourcePtr spParamSource( + new RbtParameterFileSource(Rbt::GetRbtFileName("data/scripts", config.strParamFile)) + ); + // Read the receptor parameter file + RbtParameterFileSourcePtr spRecepPrmSource( + new RbtParameterFileSource(Rbt::GetRbtFileName("data/receptors", config.strReceptorPrmFile)) + ); + cout << endl + << "DOCKING PROTOCOL:" << endl + << spParamSource->GetFileName() << endl + << spParamSource->GetTitle() << endl; + cout << endl + << "RECEPTOR:" << endl + << spRecepPrmSource->GetFileName() << endl + << spRecepPrmSource->GetTitle() << endl; + + // Create the scoring function from the SCORE section of the docking protocol prm file + // Format is: + // SECTION SCORE + // INTER RbtInterSF.prm + // INTRA RbtIntraSF.prm + // END_SECTION + // + // Notes: + // Section name must be SCORE. This is also the name of the root SF aggregate + // An aggregate is created for each parameter in the section. + // Parameter name becomes the name of the subaggregate (e.g. SCORE.INTER) + // Parameter value is the file name for the subaggregate definition + // Default directory is $RBT_ROOT/data/sf + RbtSFFactoryPtr spSFFactory(new RbtSFFactory()); // Factory class for scoring functions + RbtSFAggPtr spSF(new RbtSFAgg(_ROOT_SF)); // Root SF aggregate + spParamSource->SetSection(_ROOT_SF); + RbtStringList sfList(spParamSource->GetParameterList()); + // Loop over all parameters in the SCORE section + for (RbtStringListConstIter sfIter = sfList.begin(); sfIter != sfList.end(); sfIter++) { + // sfFile = file name for scoring function subaggregate + RbtString sfFile(Rbt::GetRbtFileName("data/sf", spParamSource->GetParameterValueAsString(*sfIter))); + RbtParameterFileSourcePtr spSFSource(new RbtParameterFileSource(sfFile)); + // Create and add the subaggregate + spSF->Add(spSFFactory->CreateAggFromFile(spSFSource, *sfIter)); + } + + // Add the RESTRAINT subaggregate scoring function from any SF definitions in the receptor prm file + spSF->Add(spSFFactory->CreateAggFromFile(spRecepPrmSource, _RESTRAINT_SF)); + + // Create the docking transform aggregate from the transform definitions in the docking prm file + RbtTransformFactoryPtr spTransformFactory(new RbtTransformFactory()); + spParamSource->SetSection(); + RbtTransformAggPtr spTransform(spTransformFactory->CreateAggFromFile(spParamSource, _ROOT_TRANSFORM)); + + // Override the TRACE levels for the scoring function and transform + // Dump details to cout + // Register the scoring function and the transform with the workspace + if (config.iTrace.has_value()) { + RbtRequestPtr spTraceReq(new RbtSFSetParamRequest("TRACE", config.iTrace.value())); + spSF->HandleRequest(spTraceReq); + spTransform->HandleRequest(spTraceReq); + } + if (config.iTrace > 0) { + cout << endl << "SCORING FUNCTION DETAILS:" << endl << *spSF << endl; + cout << endl << "SEARCH DETAILS:" << endl << *spTransform << endl; + } + spWS->SetSF(spSF); + spWS->SetTransform(spTransform); + + // DM 18 May 1999 + // Variants describing the library version, exe version, parameter file, and current directory + // Will be stored in the ligand SD files + RbtVariant vLib(Rbt::GetProduct() + " (" + Rbt::GetVersion() + ", Build" + Rbt::GetBuild() + ")"); + RbtVariant vExe(strExeName + " - " + RBT_VERSION); + RbtVariant vRecep(spRecepPrmSource->GetFileName()); + RbtVariant vPrm(spParamSource->GetFileName()); + RbtVariant vDir(Rbt::GetCurrentDirectory()); + + spRecepPrmSource->SetSection(); + // Read docking site from file and register with workspace + RbtString strASFile = spWS->GetName() + ".as"; + RbtString strInputFile = Rbt::GetRbtFileName("data/grids", strASFile); + // DM 26 Sep 2000 - ios_base::binary is invalid with IRIX CC + ifstream istr(strInputFile.c_str(), Rbt::inputMode); + // DM 14 June 2006 - bug fix to one of the longest standing rDock issues + //(the cryptic "Error reading from input stream" message, if cavity file was missing) + if (!istr) { + RbtString message = "Cavity file (" + strASFile + ") not found in current directory or $RBT_HOME"; + message += " - run rbcavity first"; + throw RbtFileReadError(_WHERE_, message); + } + RbtDockingSitePtr spDS(new RbtDockingSite(istr)); + istr.close(); + spWS->SetDockingSite(spDS); + cout << endl << "DOCKING SITE" << endl << (*spDS) << endl; + + // Prepare the SD file sink for saving the docked conformations for each ligand + // DM 3 Dec 1999 - replaced ostringstream with RbtString in determining SD file name + // SRC 2014 moved here this block to allow WRITE_ERROR TRUE + if (config.strOutputPrefix.has_value()) { + RbtMolecularFileSinkPtr spMdlFileSink(new RbtMdlFileSink(config.strOutputPrefix.value() + ".sd", RbtModelPtr()) + ); + spWS->SetSink(spMdlFileSink); + } + + RbtPRMFactory prmFactory(spRecepPrmSource, spDS); + prmFactory.SetTrace(config.iTrace.value_or(0)); + // Create the receptor model from the file names in the receptor parameter file + RbtModelPtr spReceptor = prmFactory.CreateReceptor(); + spWS->SetReceptor(spReceptor); + + // Register any solvent + RbtModelList solventList = prmFactory.CreateSolvent(); + spWS->SetSolvent(solventList); + if (spWS->hasSolvent()) { + RbtInt nSolvent = spWS->GetSolvent().size(); + cout << endl << nSolvent << " solvent molecules registered" << endl; + } else { + cout << endl << "No solvent" << endl; + } + + // SRC 2014 removed sector bOutput from here to some blocks above, for WRITEERRORS TRUE + + // Seed the random number generator + RbtRand &theRand = Rbt::GetRbtRand(); // ref to random number generator + if (config.nSeed.has_value()) { + theRand.Seed(config.nSeed.value()); + } + + // Create the filter object for controlling early termination of protocol + RbtFilterPtr spfilter; + if (config.strFilterFile.has_value()) { + spfilter = new RbtFilter(config.strFilterFile.value()); + if (config.nDockingRuns.has_value()) { + spfilter->SetMaxNRuns(config.nDockingRuns.value()); + } + } else { + spfilter = new RbtFilter(config.get_filter_string(), true); + } + if (config.iTrace.has_value()) { + RbtRequestPtr spTraceReq(new RbtSFSetParamRequest("TRACE", config.iTrace.value())); + spfilter->HandleRequest(spTraceReq); + } + + // Register the Filter with the workspace + spWS->SetFilter(spfilter); + + // MAIN LOOP OVER LIGAND RECORDS + // DM 20 Apr 1999 - add explicit bPosIonise and bNegIonise flags to MdlFileSource constructor + RbtMolecularFileSourcePtr spMdlFileSource( + new RbtMdlFileSource(config.strLigandMdlFile, config.bPosIonise, config.bNegIonise, !config.bAllH) + ); + for (RbtInt nRec = 1; spMdlFileSource->FileStatusOK(); spMdlFileSource->NextRecord(), nRec++) { + cout.setf(ios_base::left, ios_base::adjustfield); + cout << endl << "**************************************************" << endl << "RECORD #" << nRec << endl; + RbtError molStatus = spMdlFileSource->Status(); + if (!molStatus.isOK()) { + cout << endl << molStatus << endl << "************************************************" << endl; + continue; + } + + // DM 26 Jul 1999 - only read the largest segment (guaranteed to be called H) + // BGD 07 Oct 2002 - catching errors created by the ligands, + // so rbdock continues with the next one, instead of + // completely stopping + try { + spMdlFileSource->SetSegmentFilterMap(Rbt::ConvertStringToSegmentMap("H")); + + if (spMdlFileSource->isDataFieldPresent("Name")) + cout << "NAME: " << spMdlFileSource->GetDataValue("Name") << endl; + if (spMdlFileSource->isDataFieldPresent("REG_Number")) + cout << "REG_Num:" << spMdlFileSource->GetDataValue("REG_Number") << endl; + cout << std::setw(30) << "RANDOM_NUMBER_SEED:" << theRand.GetSeed() << endl; + + // Create and register the ligand model + RbtModelPtr spLigand = prmFactory.CreateLigand(spMdlFileSource); + RbtString strMolName = spLigand->GetName(); + spWS->SetLigand(spLigand); + // Update any model coords from embedded chromosomes in the ligand file + spWS->UpdateModelCoordsFromChromRecords(spMdlFileSource, config.iTrace.value_or(0)); + + // DM 18 May 1999 - store run info in model data + // Clear any previous Rbt.* data fields + spLigand->ClearAllDataFields("Rbt."); + spLigand->SetDataValue("Rbt.Library", vLib); + spLigand->SetDataValue("Rbt.Executable", vExe); + spLigand->SetDataValue("Rbt.Receptor", vRecep); + spLigand->SetDataValue("Rbt.Parameter_File", vPrm); + spLigand->SetDataValue("Rbt.Current_Directory", vDir); + + // DM 10 Dec 1999 - if in target mode, loop until target score is reached + RbtBool bTargetMet = false; + + //////////////////////////////////////////////////// + // MAIN LOOP OVER EACH SIMULATED ANNEALING RUN + // Create a history file sink, just in case it's needed by any + // of the transforms + RbtInt iRun = 1; + // need to check this here. The termination + // filter is only run once at least + // one docking run has been done. + if (config.nDockingRuns < 1) bTargetMet = true; + while (!bTargetMet) { + // Catching errors with this specific run + try { + if (config.strOutputPrefix.has_value()) { + ostringstream histr; + histr << config.strOutputPrefix.value() << "_" << strMolName << nRec << "_his_" << iRun + << ".sd"; + RbtMolecularFileSinkPtr spHistoryFileSink(new RbtMdlFileSink(histr.str(), spLigand)); + spWS->SetHistorySink(spHistoryFileSink); + } + spWS->Run(); // Dock! + RbtBool bterm = spfilter->Terminate(); + RbtBool bwrite = spfilter->Write(); + if (bterm) bTargetMet = true; + if (config.strOutputPrefix.has_value() && bwrite) { + spWS->Save(); + } + iRun++; + } catch (RbtDockingError &e) { + cout << e << endl; + } + } + // END OF MAIN LOOP OVER EACH SIMULATED ANNEALING RUN + //////////////////////////////////////////////////// + } + // END OF TRY + catch (RbtLigandError &e) { + cout << e << endl; + } + } + // END OF MAIN LOOP OVER LIGAND RECORDS + //////////////////////////////////////////////////// + cout << endl << "END OF RUN" << endl; + // if (bOutput && flexRec) { + // RbtMolecularFileSinkPtr spRecepSink(new RbtCrdFileSink(strRunName+".crd",spReceptor)); + // spRecepSink->Render(); + // } +} \ No newline at end of file diff --git a/tests/src/test_rbcavity_argparser.cpp b/tests/src/test_rbcavity_argparser.cpp index 36865f2..9904e69 100644 --- a/tests/src/test_rbcavity_argparser.cpp +++ b/tests/src/test_rbcavity_argparser.cpp @@ -1,24 +1,22 @@ #include +#include +#include #include "rbcavity/rbcavity_argparser.h" #include "rbcavity/rbcavity_config.h" - #include "test_utils.hpp" -#include -#include - - -struct TestConfig{ +struct TestConfig { std::vector input; RBCavity::RBCavityConfig expected_config; }; - TEST_CASE("rbcavity", "[argparser]") { auto config_base = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm"}; - auto config_with_list = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bList = true, .dist = 3.0f}; - auto config_with_border = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bBorder = true, .border = 5.0f}; + auto config_with_list = + RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bList = true, .dist = 3.0f}; + auto config_with_border = + RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bBorder = true, .border = 5.0f}; auto config_with_read_as = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bReadAS = true}; auto config_with_write_as = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bWriteAS = true}; auto config_with_dump_insight = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bDump = true}; @@ -26,7 +24,8 @@ TEST_CASE("rbcavity", "[argparser]") { auto config_with_site = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bSite = true}; auto config_with_moe_grid = RBCavity::RBCavityConfig{.strReceptorPrmFile = "receptor.prm", .bMOEgrid = true}; - std::vector all_settings{"-r", "receptor.prm", "-l", "3", "-b", "5", "-R", "-W", "-d", "-v", "-s", "-m"}; + std::vector all_settings{ + "-r", "receptor.prm", "-l", "3", "-b", "5", "-R", "-W", "-d", "-v", "-s", "-m"}; auto config_all_settings = RBCavity::RBCavityConfig{ .strReceptorPrmFile = "receptor.prm", .bReadAS = true, @@ -41,7 +40,7 @@ TEST_CASE("rbcavity", "[argparser]") { .dist = 3.0f, }; - std::vector configs { + std::vector configs{ TestConfig{.input = {"-r", "receptor.prm"}, .expected_config = config_base}, TestConfig{.input = {"-rreceptor.prm"}, .expected_config = config_base}, TestConfig{.input = {"-receptor=receptor.prm"}, .expected_config = config_base}, @@ -75,16 +74,14 @@ TEST_CASE("rbcavity", "[argparser]") { TestConfig{.input = {"-r", "receptor.prm", "-m"}, .expected_config = config_with_moe_grid}, TestConfig{.input = {"-r", "receptor.prm", "-dump-moe"}, .expected_config = config_with_moe_grid}, TestConfig{.input = {"-r", "receptor.prm", "--dump-moe"}, .expected_config = config_with_moe_grid}, - TestConfig{.input = all_settings, .expected_config = config_all_settings} - }; + TestConfig{.input = all_settings, .expected_config = config_all_settings}}; NullBuffer null_buffer; CoutRedirect redirect_cout(&null_buffer); CerrRedirect redirect_cerr(&null_buffer); - - for (auto test_config: configs) { - auto inputs = std::vector{"rbcavity"}; - for (auto & input: test_config.input) inputs.push_back(input.c_str()); - auto result = RBCavity::parse_args(inputs.size(), inputs.data()); - REQUIRE(result == test_config.expected_config); - } + + auto test_config = GENERATE_COPY(from_range(configs)); + auto inputs = std::vector{"rbcavity"}; + for (auto &input: test_config.input) inputs.push_back(input.c_str()); + auto result = RBCavity::parse_args(inputs.size(), inputs.data()); + REQUIRE(result == test_config.expected_config); } diff --git a/tests/src/test_rbdock_argparser.cpp b/tests/src/test_rbdock_argparser.cpp new file mode 100644 index 0000000..e9edbee --- /dev/null +++ b/tests/src/test_rbdock_argparser.cpp @@ -0,0 +1,114 @@ +#include +#include +#include + +#include "rbdock/rbdock_argparser.h" +#include "rbdock/rbdock_config.h" +#include "test_utils.hpp" + +using RBDock::RBDockConfig; + +typedef std::tuple, RBDockConfig> TestConfig; + +TEST_CASE("rbdock", "[argparser]") { + std::vector configs{ + // base config + {{"-r", "receptor.prm", "-i", "ligand.sdf"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm"}}, + // config with output + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-o", "output"}, + + RBDockConfig{ + .strLigandMdlFile = "ligand.sdf", + .strReceptorPrmFile = "receptor.prm", + .strOutputPrefix = "output", + }}, + // config with protocol + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-p", "protocol.prm"}, + + RBDockConfig{ + .strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .strParamFile = "protocol.prm"}}, + // config with filter + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-t", "filter.txt"}, + + RBDockConfig{ + .strLigandMdlFile = "ligand.sdf", + .strReceptorPrmFile = "receptor.prm", + .strFilterFile = "filter.txt", + }}, + // config with target score + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-t", "0.5"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .dTargetScore = 0.5}}, + // config with continue + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-C"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bContinue = true}}, + // config with docking runs + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-n", "5"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .nDockingRuns = 5}}, + // config with protonate + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-P"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bPosIonise = true}}, + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-ap"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bPosIonise = true}}, + // config with deprotonate + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-D"}, + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bNegIonise = true}}, + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-an"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bNegIonise = true}}, + // config with all hydrogens + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-H"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bAllH = true}}, + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-allH"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .bAllH = true}}, + // config with seed + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-s", "123"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .nSeed = 123}}, + // config with trace + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-T", "2"}, + + RBDockConfig{.strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .iTrace = 2}}, + // config with filter + {{"-r", "receptor.prm", "-i", "ligand.sdf", "-f", "filter.txt"}, + RBDockConfig{ + .strLigandMdlFile = "ligand.sdf", .strReceptorPrmFile = "receptor.prm", .strFilterFile = "filter.txt"}}, + + // config with all options + { + {"-r", "receptor.prm", "-i", "ligand.sdf", "-o", "output", "-p", "protocol.prm", "-t", "0.5", + "-C", "-n", "5", "-P", "-D", "-H", "-s", "123", "-T", "2"}, + + RBDockConfig{ + .strLigandMdlFile = "ligand.sdf", + .strReceptorPrmFile = "receptor.prm", + .strParamFile = "protocol.prm", + .strOutputPrefix = "output", + .nDockingRuns = 5, + .dTargetScore = 0.5, + .nSeed = 123, + .iTrace = 2, + .bContinue = true, + .bPosIonise = true, + .bNegIonise = true, + .bAllH = true, + }, + }, + }; + NullBuffer null_buffer; + CoutRedirect redirect_cout(&null_buffer); + CerrRedirect redirect_cerr(&null_buffer); + + auto [test_input, expected_config] = GENERATE_COPY(from_range(configs)); + + auto inputs = std::vector{"rbdock"}; + for (auto &input: test_input) inputs.push_back(input.c_str()); + std::cerr << "inputs: "; + for (auto &input: inputs) std::cerr << input << " "; + std::cerr << std::endl; + auto result = RBDock::parse_args(inputs.size(), inputs.data()); + REQUIRE(result == expected_config); +}