diff --git a/doc/colvars-code-refs.bib b/doc/colvars-code-refs.bib index 6aa885010..ff9b04273 100644 --- a/doc/colvars-code-refs.bib +++ b/doc/colvars-code-refs.bib @@ -342,4 +342,5 @@ @article{White2014 % coordNum pairlist % Custom functions (Lepton) % Scripted functions (Tcl) +% torchANN colvar component % --- END --- diff --git a/gromacs/cmake/gmxManageColvars.cmake b/gromacs/cmake/gmxManageColvars.cmake index d548d60f2..f8057f435 100644 --- a/gromacs/cmake/gmxManageColvars.cmake +++ b/gromacs/cmake/gmxManageColvars.cmake @@ -55,3 +55,15 @@ function(gmx_include_colvars_headers) target_include_directories(libgromacs PRIVATE ${PROJECT_SOURCE_DIR}/src/external/colvars) endif() endfunction() + +function(gmx_set_colvars_torch) + find_package(Torch) + if (Torch_FOUND) + set_property(TARGET colvars PROPERTY CXX_STANDARD 17) + target_compile_definitions(colvars PRIVATE -DTORCH) + target_compile_options(colvars PRIVATE ${CMAKE_CXX_FLAGS} ${TORCH_CXX_FLAGS}) + target_include_directories(colvars PRIVATE ${TORCH_INCLUDE_DIRS}) + target_link_libraries(libgromacs PRIVATE "${TORCH_LIBRARIES}") + endif() +endfunction() + diff --git a/gromacs/gromacs-2020.x.patch b/gromacs/gromacs-2020.x.patch index baef003be..af08b94b0 100644 --- a/gromacs/gromacs-2020.x.patch +++ b/gromacs/gromacs-2020.x.patch @@ -3,9 +3,9 @@ index 0911eb2a45..5530c1576a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -200,6 +200,9 @@ option(GMX_USE_OPENCL "Enable OpenCL acceleration" OFF) - + option(GMX_INSTALL_LEGACY_API "Install legacy headers" OFF) - + +include(gmxManageColvars) +include(gmxManageLepton) + @@ -13,13 +13,13 @@ index 0911eb2a45..5530c1576a 100644 set(REQUIRED_CUDA_VERSION 9.0) set(REQUIRED_CUDA_COMPUTE_CAPABILITY 3.0) diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt -index 9249a7a08f..18b45ff01f 100644 +index 9249a7a08f..ab5a927fb7 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -137,6 +137,12 @@ if (WIN32) endif() list(APPEND libgromacs_object_library_dependencies thread_mpi) - + +# Add Colvars and Lepton targets, embed their object code in libgromacs +gmx_manage_colvars() +gmx_manage_lepton() @@ -29,11 +29,13 @@ index 9249a7a08f..18b45ff01f 100644 configure_file(version.h.cmakein version.h) if(GMX_INSTALL_LEGACY_API) install(FILES -@@ -195,6 +201,8 @@ else() +@@ -195,6 +201,10 @@ else() add_library(libgromacs ${LIBGROMACS_SOURCES}) endif() - + +gmx_include_colvars_headers() ++ ++gmx_set_colvars_torch() + # Add these contents first because linking their tests can take a lot # of time, so we want lots of parallel work still available after @@ -58,8 +60,8 @@ index 9c6cfe4213..ec2e64f61f 100644 +/* COLVARS : add a value of 1000 for Colvars support and + * prevent regular GROMACS to read colvars .cpt files */ +static const int cpt_version = cptv_Count - 1 + 1000; - - + + const char* est_names[estNR] = { "FE-lambda", @@ -1178,6 +1181,15 @@ static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderC { @@ -75,12 +77,12 @@ index 9c6cfe4213..ec2e64f61f 100644 + } + } - + static int do_cpt_footer(XDR* xd, int file_version) @@ -1909,6 +1921,35 @@ static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDst return 0; } - + +/* This function stores the last whole configuration of the colvars atoms in the .cpt file */ +static int do_cpt_colvars(XDR* xd, gmx_bool bRead, int ecolvars, colvarshistory_t* colvarsstate, FILE* list) +{ @@ -116,7 +118,7 @@ index 9c6cfe4213..ec2e64f61f 100644 @@ -2330,6 +2371,10 @@ void write_checkpoint(const char* fn, swaphistory_t* swaphist = observablesHistory->swapHistory.get(); int eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO); - + + /* COLVARS */ + colvarshistory_t* colvarshist = observablesHistory->colvarsHistory.get(); + int ecolvars = (colvarshist ? colvarshist->n_atoms : 0); @@ -145,7 +147,7 @@ index 9c6cfe4213..ec2e64f61f 100644 @@ -2802,6 +2849,17 @@ static void read_checkpoint(const char* fn, cp_error(); } - + + if (headerContents->ecolvars != 0 && observablesHistory->colvarsHistory == nullptr) + { + observablesHistory->colvarsHistory = std::make_unique(colvarshistory_t{}); @@ -163,7 +165,7 @@ index 9c6cfe4213..ec2e64f61f 100644 @@ -2957,6 +3015,13 @@ static CheckpointHeaderContents read_checkpoint_data(t_fileio* cp_error(); } - + + colvarshistory_t colvarshist = {}; + ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, headerContents.ecolvars, &colvarshist, nullptr); + if (ret) @@ -172,12 +174,12 @@ index 9c6cfe4213..ec2e64f61f 100644 + } + ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version); - + if (ret) @@ -3065,6 +3130,12 @@ void list_checkpoint(const char* fn, FILE* out) ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out); } - + + if (ret == 0) + { + colvarshistory_t colvarshist = {}; @@ -198,7 +200,7 @@ index fb8f7268be..6feb181b30 100644 + //! Colvars + int ecolvars; }; - + /* Write a checkpoint to .cpt diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index f2532f3dfe..c761958854 100644 @@ -210,7 +212,7 @@ index f2532f3dfe..c761958854 100644 bEner_[F_ORIRESDEV] = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0); - bEner_[F_COM_PULL] = ((ir->bPull && pull_have_potential(pull_work)) || ir->bRot); + bEner_[F_COM_PULL] = ((ir->bPull && pull_have_potential(pull_work)) || ir->bRot || ir->bColvars); - + MdModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest; mdModulesNotifier.notifier_.notify(&mdModulesAddOutputToDensityFittingFieldRequest); diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp @@ -220,7 +222,7 @@ index f2528d78b4..7f91cac83d 100644 @@ -114,6 +114,8 @@ #include "gromacs/utility/strconvert.h" #include "gromacs/utility/sysinfo.h" - + +#include "colvarproxy_gromacs.h" + using gmx::AtomLocality; @@ -242,7 +244,7 @@ index f2528d78b4..7f91cac83d 100644 + gmx::ForceProviderInput forceProviderInput(x, *mdatoms, t, box, *cr); gmx::ForceProviderOutput forceProviderOutput(forceWithVirial, enerd); - + diff --git a/src/gromacs/mdrun/legacymdrunoptions.h b/src/gromacs/mdrun/legacymdrunoptions.h index 796e479490..8b20073b3b 100644 --- a/src/gromacs/mdrun/legacymdrunoptions.h @@ -255,7 +257,7 @@ index 796e479490..8b20073b3b 100644 + { efXVG, "-swap", "swapions", ffOPTWR }, + { efDAT, "-colvars", "colvars", ffOPTRDMULT }, /* COLVARS */ + { efDAT, "-colvars_restart", "colvars", ffOPTRD }, /* COLVARS */}}; - + //! Print a warning if any force is larger than this (in kJ/mol nm). real pforce = -1; diff --git a/src/gromacs/mdrun/replicaexchange.cpp b/src/gromacs/mdrun/replicaexchange.cpp @@ -268,10 +270,10 @@ index 9ff4b3817d..eb31f1fa89 100644 exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms); + exchange_rvecs(ms, b, state->xa_old_whole_colvars, state->n_colvars_atoms); } - + static void copy_state_serial(const t_state* src, t_state* dest) diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp -index c2b3c088d7..bb38d44ab2 100644 +index c2b3c088d7..7a5f2b05a4 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -115,6 +115,7 @@ @@ -285,7 +287,7 @@ index c2b3c088d7..bb38d44ab2 100644 @@ -156,6 +157,8 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringutil.h" - + +#include "colvarproxy_gromacs.h" + #include "isimulator.h" @@ -294,7 +296,7 @@ index c2b3c088d7..bb38d44ab2 100644 @@ -1536,6 +1539,51 @@ int Mdrunner::mdrunner() MASTER(cr) ? globalState->x.rvec_array() : nullptr, filenames.size(), filenames.data(), oenv, mdrunOptions.imdOptions, startingBehavior); - + + /* COLVARS */ + if (opt2bSet("-colvars",filenames.size(), filenames.data())) + { @@ -346,7 +348,7 @@ index c2b3c088d7..bb38d44ab2 100644 @@ -1654,6 +1702,16 @@ int Mdrunner::mdrunner() free_membed(membed); } - + + /* COLVARS */ + if (inputrec->bColvars) + { @@ -387,20 +389,20 @@ index 0000000000..ea69e03419 + +#endif diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 266670f3ef..f1d287d35d 100644 +index 266670f3ef..86e4aa35d3 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -53,6 +53,8 @@ struct gmx_enfrot; struct gmx_enfrotgrp; struct pull_params_t; - + +class colvarproxy_gromacs; + namespace gmx { class Awh; @@ -587,6 +589,10 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding) - + //! KVT for storing simulation parameters that are not part of the mdp file. std::unique_ptr internalParameters; + @@ -408,7 +410,7 @@ index 266670f3ef..f1d287d35d 100644 + gmx_bool bColvars = false; /* Do we do colvars calculations ? */ + colvarproxy_gromacs *colvars_proxy = nullptr; /* The object for the colvars calculations */ }; - + int ir_optimal_nstcalcenergy(const t_inputrec* ir); diff --git a/src/gromacs/mdtypes/observableshistory.cpp b/src/gromacs/mdtypes/observableshistory.cpp index 0b5983a59c..57d851645a 100644 @@ -419,7 +421,7 @@ index 0b5983a59c..57d851645a 100644 #include "gromacs/mdtypes/pullhistory.h" #include "gromacs/mdtypes/swaphistory.h" +#include "gromacs/mdtypes/colvarshistory.h" - + ObservablesHistory::ObservablesHistory() = default; ObservablesHistory::~ObservablesHistory() = default; diff --git a/src/gromacs/mdtypes/observableshistory.h b/src/gromacs/mdtypes/observableshistory.h @@ -431,18 +433,18 @@ index d2ba1d820f..a5747139d7 100644 struct edsamhistory_t; struct swaphistory_t; +struct colvarshistory_t; - + /*! \libinternal \brief Observables history, for writing/reading to/from checkpoint file */ @@ -76,6 +77,9 @@ struct ObservablesHistory //! Ion/water position swapping history std::unique_ptr swapHistory; - + + //! Colvars + std::unique_ptr colvarsHistory; + ObservablesHistory(); - + ~ObservablesHistory(); diff --git a/src/gromacs/mdtypes/state.cpp b/src/gromacs/mdtypes/state.cpp index a949d589a3..957f7b7139 100644 @@ -456,7 +458,7 @@ index a949d589a3..957f7b7139 100644 + ddp_count_cg_gl(0), + xa_old_whole_colvars(nullptr), + n_colvars_atoms(0) - + { // It would be nicer to initialize these with {} or {{0}} in the diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h @@ -465,14 +467,14 @@ index a54bff29bb..8619c7935a 100644 +++ b/src/gromacs/mdtypes/state.h @@ -257,6 +257,10 @@ public: std::vector cg_gl; //!< The global cg number of the local cgs - + std::vector pull_com_prev_step; //!< The COM of the previous step of each pull group + + int n_colvars_atoms; //!< number of colvars atoms + rvec* xa_old_whole_colvars; //!< last whole positions of colvars atoms + }; - + #ifndef DOXYGEN diff --git a/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml b/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml index c2973bb1af..cb4d1da254 100644 @@ -496,6 +498,6 @@ index c2973bb1af..cb4d1da254 100644 + Generic data file + -colvars_restart [<.dat>] (colvars.dat) (Opt.) + Generic data file - + Options to specify output files: - + diff --git a/gromacs/gromacs-2021.x.patch b/gromacs/gromacs-2021.x.patch index 3f8bdf9f9..7b2f82710 100644 --- a/gromacs/gromacs-2021.x.patch +++ b/gromacs/gromacs-2021.x.patch @@ -3,9 +3,9 @@ index 3715ce1656..65134d0568 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -572,6 +572,9 @@ include(gmxManageTNG) - + include(gmxManageLmfit) - + +include(gmxManageColvars) +include(gmxManageLepton) + @@ -13,13 +13,13 @@ index 3715ce1656..65134d0568 100644 set(REQUIRED_CUDA_VERSION 9.0) set(REQUIRED_CUDA_COMPUTE_CAPABILITY 3.0) diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt -index a4430e9dd6..385b688501 100644 +index a4430e9dd6..74dc70817e 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -142,6 +142,12 @@ if (WIN32) endif() list(APPEND libgromacs_object_library_dependencies thread_mpi) - + +# Add Colvars and Lepton targets, embed their object code in libgromacs +gmx_manage_colvars() +gmx_manage_lepton() @@ -29,11 +29,13 @@ index a4430e9dd6..385b688501 100644 configure_file(version.h.cmakein version.h) if(GMX_INSTALL_LEGACY_API) install(FILES -@@ -189,6 +195,8 @@ else() +@@ -189,6 +195,10 @@ else() add_library(libgromacs ${LIBGROMACS_SOURCES}) endif() - + +gmx_include_colvars_headers() ++ ++gmx_set_colvars_torch() + # Add these contents first because linking their tests can take a lot # of time, so we want lots of parallel work still available after @@ -58,8 +60,8 @@ index fadc3d283e..8a6ce56a6f 100644 +/* COLVARS : add a value of 1000 for Colvars support and + * prevent regular GROMACS to read colvars .cpt files */ +static const int cpt_version = cptv_Count - 1 + 1000; - - + + const char* est_names[estNR] = { "FE-lambda", @@ -1238,6 +1241,15 @@ static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderC { @@ -75,12 +77,12 @@ index fadc3d283e..8a6ce56a6f 100644 + } + } - + static int do_cpt_footer(XDR* xd, int file_version) @@ -1964,6 +1976,35 @@ static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDst return 0; } - + +/* This function stores the last whole configuration of the colvars atoms in the .cpt file */ +static int do_cpt_colvars(XDR* xd, gmx_bool bRead, int ecolvars, colvarshistory_t* colvarsstate, FILE* list) +{ @@ -124,7 +126,7 @@ index fadc3d283e..8a6ce56a6f 100644 @@ -2713,6 +2755,17 @@ static void read_checkpoint(const char* fn, cp_error(); } - + + if (headerContents->ecolvars != 0 && observablesHistory->colvarsHistory == nullptr) + { + observablesHistory->colvarsHistory = std::make_unique(colvarshistory_t{}); @@ -142,7 +144,7 @@ index fadc3d283e..8a6ce56a6f 100644 @@ -2879,6 +2932,13 @@ static CheckpointHeaderContents read_checkpoint_data(t_fileio* cp_error(); } - + + colvarshistory_t colvarshist = {}; + ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, headerContents.ecolvars, &colvarshist, nullptr); + if (ret) @@ -151,12 +153,12 @@ index fadc3d283e..8a6ce56a6f 100644 + } + ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version); - + if (ret) @@ -3000,6 +3060,12 @@ void list_checkpoint(const char* fn, FILE* out) ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out); } - + + if (ret == 0) + { + colvarshistory_t colvarshist = {}; @@ -189,7 +191,7 @@ index 0b72fc4e0c..60666542b4 100644 bEner_[F_ORIRESDEV] = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0); - bEner_[F_COM_PULL] = ((ir->bPull && pull_have_potential(*pull_work)) || ir->bRot); + bEner_[F_COM_PULL] = ((ir->bPull && pull_have_potential(*pull_work)) || ir->bRot || ir->bColvars); - + MdModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest; mdModulesNotifier.simulationSetupNotifications_.notify(&mdModulesAddOutputToDensityFittingFieldRequest); diff --git a/src/gromacs/mdlib/mdoutf.cpp b/src/gromacs/mdlib/mdoutf.cpp @@ -207,7 +209,7 @@ index 7e06e4fc9e..38f5604ebd 100644 @@ -357,6 +358,10 @@ static void write_checkpoint(const char* fn, swaphistory_t* swaphist = observablesHistory->swapHistory.get(); int eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO); - + + /* COLVARS */ + colvarshistory_t* colvarshist = observablesHistory->colvarsHistory.get(); + int ecolvars = (colvarshist ? colvarshist->n_atoms : 0); @@ -230,11 +232,11 @@ index 2571b0d216..4572e2baea 100644 @@ -123,6 +123,8 @@ #include "gromacs/utility/strconvert.h" #include "gromacs/utility/sysinfo.h" - + +#include "colvarproxy_gromacs.h" + #include "gpuforcereduction.h" - + using gmx::ArrayRef; @@ -616,6 +618,16 @@ static void computeSpecialForces(FILE* fplog, */ @@ -252,7 +254,7 @@ index 2571b0d216..4572e2baea 100644 + gmx::ForceProviderInput forceProviderInput(x, *mdatoms, t, box, *cr); gmx::ForceProviderOutput forceProviderOutput(forceWithVirialMtsLevel0, enerd); - + diff --git a/src/gromacs/mdrun/legacymdrunoptions.h b/src/gromacs/mdrun/legacymdrunoptions.h index 474f6f0396..bb94199a30 100644 --- a/src/gromacs/mdrun/legacymdrunoptions.h @@ -265,7 +267,7 @@ index 474f6f0396..bb94199a30 100644 + { efXVG, "-swap", "swapions", ffOPTWR }, + { efDAT, "-colvars", "colvars", ffOPTRDMULT }, /* COLVARS */ + { efDAT, "-colvars_restart", "colvars", ffOPTRD }, /* COLVARS */}}; - + //! Print a warning if any force is larger than this (in kJ/mol nm). real pforce = -1; diff --git a/src/gromacs/mdrun/replicaexchange.cpp b/src/gromacs/mdrun/replicaexchange.cpp @@ -278,10 +280,10 @@ index c40161d9ef..490db3f10f 100644 exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms); + exchange_rvecs(ms, b, state->xa_old_whole_colvars, state->n_colvars_atoms); } - + static void copy_state_serial(const t_state* src, t_state* dest) diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp -index 232d994e1a..8937b83296 100644 +index 232d994e1a..cdf06fe041 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -125,6 +125,7 @@ @@ -295,7 +297,7 @@ index 232d994e1a..8937b83296 100644 @@ -167,6 +168,8 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringutil.h" - + +#include "colvarproxy_gromacs.h" + #include "isimulator.h" @@ -304,7 +306,7 @@ index 232d994e1a..8937b83296 100644 @@ -1691,6 +1694,51 @@ int Mdrunner::mdrunner() MASTER(cr) ? globalState->x.rvec_array() : nullptr, filenames.size(), filenames.data(), oenv, mdrunOptions.imdOptions, startingBehavior); - + + /* COLVARS */ + if (opt2bSet("-colvars",filenames.size(), filenames.data())) + { @@ -356,7 +358,7 @@ index 232d994e1a..8937b83296 100644 @@ -1839,6 +1887,16 @@ int Mdrunner::mdrunner() } releaseDevice(deviceInfo); - + + /* COLVARS */ + if (inputrec->bColvars) + { @@ -397,20 +399,20 @@ index 0000000000..6605e6fce2 + +#endif diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 6e4ee727ab..042acdc01b 100644 +index 6e4ee727ab..f827185107 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -55,6 +55,8 @@ struct gmx_enfrot; struct gmx_enfrotgrp; struct pull_params_t; - + +class colvarproxy_gromacs; + namespace gmx { class Awh; @@ -570,6 +572,10 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding) - + //! KVT for storing simulation parameters that are not part of the mdp file. std::unique_ptr internalParameters; + @@ -418,7 +420,7 @@ index 6e4ee727ab..042acdc01b 100644 + gmx_bool bColvars = false; /* Do we do colvars calculations ? */ + colvarproxy_gromacs *colvars_proxy = nullptr; /* The object for the colvars calculations */ }; - + int ir_optimal_nstcalcenergy(const t_inputrec* ir); diff --git a/src/gromacs/mdtypes/observableshistory.cpp b/src/gromacs/mdtypes/observableshistory.cpp index 0b5983a59c..57d851645a 100644 @@ -429,7 +431,7 @@ index 0b5983a59c..57d851645a 100644 #include "gromacs/mdtypes/pullhistory.h" #include "gromacs/mdtypes/swaphistory.h" +#include "gromacs/mdtypes/colvarshistory.h" - + ObservablesHistory::ObservablesHistory() = default; ObservablesHistory::~ObservablesHistory() = default; diff --git a/src/gromacs/mdtypes/observableshistory.h b/src/gromacs/mdtypes/observableshistory.h @@ -441,18 +443,18 @@ index d2ba1d820f..a5747139d7 100644 struct edsamhistory_t; struct swaphistory_t; +struct colvarshistory_t; - + /*! \libinternal \brief Observables history, for writing/reading to/from checkpoint file */ @@ -76,6 +77,9 @@ struct ObservablesHistory //! Ion/water position swapping history std::unique_ptr swapHistory; - + + //! Colvars + std::unique_ptr colvarsHistory; + ObservablesHistory(); - + ~ObservablesHistory(); diff --git a/src/gromacs/mdtypes/state.cpp b/src/gromacs/mdtypes/state.cpp index 0f36009513..b42ba3caf7 100644 @@ -466,7 +468,7 @@ index 0f36009513..b42ba3caf7 100644 + ddp_count_cg_gl(0), + xa_old_whole_colvars(nullptr), + n_colvars_atoms(0) - + { // It would be nicer to initialize these with {} or {{0}} in the diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h @@ -475,13 +477,13 @@ index e38f3f7dbc..06a1cd8484 100644 +++ b/src/gromacs/mdtypes/state.h @@ -269,6 +269,9 @@ public: std::vector cg_gl; //!< The global cg number of the local cgs - + std::vector pull_com_prev_step; //!< The COM of the previous step of each pull group + + int n_colvars_atoms; //!< number of colvars atoms + rvec* xa_old_whole_colvars; //!< last whole positions of colvars atoms }; - + #ifndef DOXYGEN diff --git a/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml b/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml index c2973bb1af..cb4d1da254 100644 @@ -505,5 +507,6 @@ index c2973bb1af..cb4d1da254 100644 + Generic data file + -colvars_restart [<.dat>] (colvars.dat) (Opt.) + Generic data file - + Options to specify output files: + diff --git a/gromacs/gromacs-2022.x.patch b/gromacs/gromacs-2022.x.patch index c4cf77ab2..9f6e3c394 100644 --- a/gromacs/gromacs-2022.x.patch +++ b/gromacs/gromacs-2022.x.patch @@ -1,11 +1,11 @@ diff --git a/CMakeLists.txt b/CMakeLists.txt -index cd748c9..31e3b62 100644 +index cd748c9..dffee16 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -588,6 +588,10 @@ include(gmxManageLmfit) - + include(gmxManageMuparser) - + +include(gmxManageColvars) + +include(gmxManageLepton) @@ -14,13 +14,13 @@ index cd748c9..31e3b62 100644 # Process SIMD instruction settings ################################################## diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt -index 99cc804..40fbe29 100644 +index 99cc804..0879cd4 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -139,6 +139,11 @@ if (WIN32) endif() list(APPEND libgromacs_object_library_dependencies thread_mpi) - + +gmx_manage_colvars() +gmx_manage_lepton() +list(APPEND libgromacs_object_library_dependencies colvars) @@ -29,17 +29,19 @@ index 99cc804..40fbe29 100644 if(GMX_INSTALL_LEGACY_API) install(FILES analysisdata.h -@@ -184,6 +189,8 @@ else() +@@ -184,6 +189,10 @@ else() add_library(libgromacs ${LIBGROMACS_SOURCES}) endif() - + +gmx_include_colvars_headers() ++ ++gmx_set_colvars_torch() + if (TARGET Heffte::Heffte) target_link_libraries(libgromacs PRIVATE Heffte::Heffte) endif() diff --git a/src/gromacs/applied_forces/CMakeLists.txt b/src/gromacs/applied_forces/CMakeLists.txt -index 34114cd..b469541 100644 +index 34114cd..8e4b3cf 100644 --- a/src/gromacs/applied_forces/CMakeLists.txt +++ b/src/gromacs/applied_forces/CMakeLists.txt @@ -73,6 +73,7 @@ gmx_add_libgromacs_sources( @@ -47,11 +49,11 @@ index 34114cd..b469541 100644 add_subdirectory(densityfitting) add_subdirectory(qmmm) +add_subdirectory(colvars) - + if (BUILD_TESTING) add_subdirectory(tests) diff --git a/src/gromacs/fileio/checkpoint.cpp b/src/gromacs/fileio/checkpoint.cpp -index f6764a1..279953d 100644 +index f6764a1..7dfc52c 100644 --- a/src/gromacs/fileio/checkpoint.cpp +++ b/src/gromacs/fileio/checkpoint.cpp @@ -69,6 +69,7 @@ @@ -62,7 +64,7 @@ index f6764a1..279953d 100644 #include "gromacs/modularsimulator/modularsimulator.h" #include "gromacs/trajectory/trajectoryframe.h" #include "gromacs/utility/arrayref.h" -@@ -1287,6 +1293,14 @@ static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderC +@@ -1287,6 +1288,14 @@ static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderC { contents->isModularSimulatorCheckpoint = false; } @@ -75,12 +77,12 @@ index f6764a1..279953d 100644 + contents->ecolvars = 0; + } } - + static int do_cpt_footer(XDR* xd, CheckPointVersion file_version) -@@ -2035,6 +2049,36 @@ static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDst +@@ -2035,6 +2044,36 @@ static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDst return 0; } - + +/* This function stores the last whole configuration of the colvars atoms in the .cpt file */ +static int do_cpt_colvars(XDR* xd, gmx_bool bRead, int ecolvars, colvarshistory_t* colvarsstate, FILE* list) +{ @@ -114,7 +116,7 @@ index f6764a1..279953d 100644 static int do_cpt_correlation_grid(XDR* xd, gmx_bool bRead, gmx_unused int fflags, -@@ -2432,6 +2476,7 @@ void write_checkpoint_data(t_fileio* fp, +@@ -2432,6 +2471,7 @@ void write_checkpoint_data(t_fileio* fp, observablesHistory->swapHistory.get(), nullptr) < 0) @@ -122,10 +124,10 @@ index f6764a1..279953d 100644 || (do_cpt_files(gmx_fio_getxdr(fp), FALSE, outputfiles, nullptr, headerContents.file_version) < 0)) { gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?"); -@@ -2825,6 +2870,17 @@ static void read_checkpoint(const char* fn, +@@ -2825,6 +2865,17 @@ static void read_checkpoint(const char* fn, cp_error(); } - + + if (headerContents->ecolvars != 0 && observablesHistory->colvarsHistory == nullptr) + { + observablesHistory->colvarsHistory = std::make_unique(colvarshistory_t{}); @@ -140,10 +142,10 @@ index f6764a1..279953d 100644 std::vector outputfiles; ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version); if (ret) -@@ -3000,6 +3056,13 @@ static CheckpointHeaderContents read_checkpoint_data(t_fileio* +@@ -3000,6 +3051,13 @@ static CheckpointHeaderContents read_checkpoint_data(t_fileio* cp_error(); } - + + colvarshistory_t colvarshist = {}; + ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, headerContents.ecolvars, &colvarshist, nullptr); + if (ret) @@ -152,12 +154,12 @@ index f6764a1..279953d 100644 + } + ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version); - + if (ret) -@@ -3120,6 +3183,12 @@ void list_checkpoint(const char* fn, FILE* out) +@@ -3120,6 +3178,12 @@ void list_checkpoint(const char* fn, FILE* out) ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out); } - + + if (ret == 0) + { + colvarshistory_t colvarshist = {}; @@ -199,7 +201,7 @@ index 6c0fca5..e06e511 100644 bEner_[F_ORIRESDEV] = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0); - bEner_[F_COM_PULL] = ((inputrec.bPull && pull_have_potential(*pull_work)) || inputrec.bRot); + bEner_[F_COM_PULL] = ((inputrec.bPull && pull_have_potential(*pull_work)) || inputrec.bRot || inputrec.bColvars); - + // Check MDModules for any energy output MDModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest; diff --git a/src/gromacs/mdlib/mdoutf.cpp b/src/gromacs/mdlib/mdoutf.cpp @@ -217,7 +219,7 @@ index 0a8653a..b75c8b1 100644 @@ -353,6 +354,10 @@ static void write_checkpoint(const char* fn, swaphistory_t* swaphist = observablesHistory->swapHistory.get(); SwapType eSwapCoords = (swaphist ? swaphist->eSwapCoords : SwapType::No); - + + /* COLVARS */ + colvarshistory_t* colvarshist = observablesHistory->colvarsHistory.get(); + int ecolvars = (colvarshist ? colvarshist->n_atoms : 0); @@ -234,17 +236,17 @@ index 0a8653a..b75c8b1 100644 std::strcpy(headerContents.version, gmx_version()); std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath()); diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp -index 5551227..a104775 100644 +index 5551227..944eed4 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -122,6 +122,8 @@ #include "gromacs/utility/stringutil.h" #include "gromacs/utility/sysinfo.h" - + +#include "gromacs/applied_forces/colvars/colvarproxy_gromacs.h" + #include "gpuforcereduction.h" - + using gmx::ArrayRef; @@ -648,6 +650,16 @@ static void computeSpecialForces(FILE* fplog, */ @@ -275,7 +277,7 @@ index 20d793e..b37d822 100644 + { efXVG, "-swap", "swapions", ffOPTWR }, + { efDAT, "-colvars", "colvars", ffOPTRDMULT }, /* COLVARS */ + { efDAT, "-colvars_restart", "colvars", ffOPTRD }, /* COLVARS */}}; - + //! Print a warning if any force is larger than this (in kJ/mol nm). real pforce = -1; diff --git a/src/gromacs/mdrun/replicaexchange.cpp b/src/gromacs/mdrun/replicaexchange.cpp @@ -288,10 +290,10 @@ index 9fa553c..9c5be82 100644 exchange_rvecs(ms, b, state->v.rvec_array(), state->natoms); + exchange_rvecs(ms, b, state->xa_old_whole_colvars, state->n_colvars_atoms); } - + static void copy_state_serial(const t_state* src, t_state* dest) diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp -index ea7f402..d3ffd59 100644 +index ea7f402..3bd9641 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -126,6 +126,7 @@ @@ -305,7 +307,7 @@ index ea7f402..d3ffd59 100644 @@ -170,6 +171,8 @@ #include "gromacs/utility/stringutil.h" #include "gromacs/utility/mpiinfo.h" - + +#include "gromacs/applied_forces/colvars/colvarproxy_gromacs.h" + #include "isimulator.h" @@ -314,7 +316,7 @@ index ea7f402..d3ffd59 100644 @@ -2055,6 +2058,52 @@ int Mdrunner::mdrunner() mdrunOptions.imdOptions, startingBehavior); - + + /* COLVARS */ + if (opt2bSet("-colvars",filenames.size(), filenames.data())) + { @@ -367,7 +369,7 @@ index ea7f402..d3ffd59 100644 @@ -2225,6 +2274,16 @@ int Mdrunner::mdrunner() releaseDevice(deviceInfo); } - + + /* COLVARS */ + if (inputrec->bColvars) + { @@ -383,7 +385,7 @@ index ea7f402..d3ffd59 100644 walltime_accounting_destroy(walltime_accounting); diff --git a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h new file mode 100644 -index 0000000000..6605e6fce2 +index 0000000..6605e6f --- /dev/null +++ b/src/gromacs/mdtypes/colvarshistory.h @@ -0,0 +1,20 @@ @@ -408,20 +410,20 @@ index 0000000000..6605e6fce2 + +#endif diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 922a22e..9d981ff 100644 +index 922a22e..ec34a6e 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -51,6 +51,8 @@ struct gmx_enfrot; struct gmx_enfrotgrp; struct pull_params_t; - + +class colvarproxy_gromacs; + namespace gmx { class Awh; @@ -577,6 +579,10 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding) - + //! KVT for storing simulation parameters that are not part of the mdp file. std::unique_ptr internalParameters; + @@ -429,7 +431,7 @@ index 922a22e..9d981ff 100644 + bool bColvars = false; /* Do we do colvars calculations ? */ + colvarproxy_gromacs *colvars_proxy = nullptr; /* The object for the colvars calculations */ }; - + int ir_optimal_nstcalcenergy(const t_inputrec* ir); diff --git a/src/gromacs/mdtypes/observableshistory.cpp b/src/gromacs/mdtypes/observableshistory.cpp index 12230ee..13f5df1 100644 @@ -440,7 +442,7 @@ index 12230ee..13f5df1 100644 #include "gromacs/mdtypes/pullhistory.h" #include "gromacs/mdtypes/swaphistory.h" +#include "gromacs/mdtypes/colvarshistory.h" - + ObservablesHistory::ObservablesHistory() = default; ObservablesHistory::~ObservablesHistory() = default; diff --git a/src/gromacs/mdtypes/observableshistory.h b/src/gromacs/mdtypes/observableshistory.h @@ -452,18 +454,18 @@ index 6ed0e3f..172ead0 100644 struct edsamhistory_t; struct swaphistory_t; +struct colvarshistory_t; - + /*! \libinternal \brief Observables history, for writing/reading to/from checkpoint file */ @@ -75,6 +76,9 @@ struct ObservablesHistory //! Ion/water position swapping history std::unique_ptr swapHistory; - + + //! Colvars + std::unique_ptr colvarsHistory; + ObservablesHistory(); - + ~ObservablesHistory(); diff --git a/src/gromacs/mdtypes/state.cpp b/src/gromacs/mdtypes/state.cpp index 99b315e..d49a386 100644 @@ -477,7 +479,7 @@ index 99b315e..d49a386 100644 + ddp_count_cg_gl(0), + xa_old_whole_colvars(nullptr), + n_colvars_atoms(0) - + { // It would be nicer to initialize these with {} or {{0}} in the diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h @@ -486,13 +488,13 @@ index 579e343..725e841 100644 +++ b/src/gromacs/mdtypes/state.h @@ -282,6 +282,9 @@ public: std::vector cg_gl; //!< The global cg number of the local cgs - + std::vector pull_com_prev_step; //!< The COM of the previous step of each pull group + + int n_colvars_atoms; //!< number of colvars atoms + rvec* xa_old_whole_colvars; //!< last whole positions of colvars atoms }; - + #ifndef DOXYGEN diff --git a/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml b/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml index c2973bb..cb4d1da 100644 @@ -516,6 +518,6 @@ index c2973bb..cb4d1da 100644 + Generic data file + -colvars_restart [<.dat>] (colvars.dat) (Opt.) + Generic data file - + Options to specify output files: - + diff --git a/gromacs/gromacs-2023.x.patch b/gromacs/gromacs-2023.x.patch index ef5a373b6..b70ac7367 100644 --- a/gromacs/gromacs-2023.x.patch +++ b/gromacs/gromacs-2023.x.patch @@ -38,7 +38,7 @@ index 8ddc9c1..21b3221 100644 int tcouple_min_integration_steps(TemperatureCoupling etc); diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt -index 114bfd9..bb03eb4 100644 +index 114bfd9..acbdc9c 100644 --- a/src/gromacs/CMakeLists.txt +++ b/src/gromacs/CMakeLists.txt @@ -146,6 +146,11 @@ if (WIN32) @@ -53,11 +53,12 @@ index 114bfd9..bb03eb4 100644 # This code is here instead of utility/CMakeLists.txt, because CMake # custom commands and source file properties can only be set in the directory # that contains the target that uses them. -@@ -192,6 +197,8 @@ else() +@@ -192,6 +197,9 @@ else() add_library(libgromacs ${LIBGROMACS_SOURCES}) endif() +gmx_include_colvars_headers() ++gmx_set_colvars_torch() + if (TARGET Heffte::Heffte) target_link_libraries(libgromacs PRIVATE Heffte::Heffte) @@ -398,7 +399,7 @@ index acadeab..faf31b5 100644 /* Does what it says */ print_date_and_time(fplog, cr->nodeid, "Finished mdrun", gmx_gettime()); walltime_accounting_destroy(walltime_accounting); -diff --git b/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h +diff --git a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h new file mode 100644 index 0000000..6605e6f --- /dev/null diff --git a/lammps/COLVARS.cmake b/lammps/COLVARS.cmake new file mode 100644 index 000000000..ddce71d11 --- /dev/null +++ b/lammps/COLVARS.cmake @@ -0,0 +1,48 @@ +set(COLVARS_SOURCE_DIR ${LAMMPS_LIB_SOURCE_DIR}/colvars) + +file(GLOB COLVARS_SOURCES CONFIGURE_DEPENDS ${COLVARS_SOURCE_DIR}/[^.]*.cpp) + +option(COLVARS_DEBUG "Enable debugging messages for Colvars (quite verbose)" OFF) + +option(COLVARS_LEPTON "Use the Lepton library for custom expressions" ON) + +option(COLVARS_TORCH "Enable torchann colvar components" OFF) + +if(COLVARS_LEPTON) + if(NOT LEPTON_SOURCE_DIR) + include(Packages/LEPTON) + endif() +endif() + + +add_library(colvars STATIC ${COLVARS_SOURCES}) +target_compile_definitions(colvars PRIVATE -DCOLVARS_LAMMPS) +separate_arguments(CMAKE_TUNE_FLAGS) +foreach(_FLAG ${CMAKE_TUNE_FLAGS}) + target_compile_options(colvars PRIVATE ${_FLAG}) +endforeach() +set_target_properties(colvars PROPERTIES OUTPUT_NAME lammps_colvars${LAMMPS_MACHINE}) +target_include_directories(colvars PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/colvars) +# The line below is needed to locate math_eigen_impl.h +target_include_directories(colvars PRIVATE ${LAMMPS_SOURCE_DIR}) +target_link_libraries(lammps PRIVATE colvars) + +if(COLVARS_DEBUG) + # Need to export the define publicly to be valid in interface code + target_compile_definitions(colvars PUBLIC -DCOLVARS_DEBUG) +endif() + +if(COLVARS_LEPTON) + target_compile_definitions(colvars PRIVATE -DLEPTON) + target_link_libraries(colvars PUBLIC lepton) +endif() + +if(COLVARS_TORCH) + find_package(Torch REQUIRED) + set_property(TARGET colvars PROPERTY CXX_STANDARD 17) + target_compile_definitions(colvars PRIVATE -DTORCH) + target_compile_options(colvars PRIVATE ${CMAKE_CXX_FLAGS} ${TORCH_CXX_FLAGS}) + target_include_directories(colvars PRIVATE ${TORCH_INCLUDE_DIRS}) + target_link_libraries(lmp PRIVATE ${TORCH_LIBRARIES}) +endif() + diff --git a/lammps/lib/colvars/Makefile.common b/lammps/lib/colvars/Makefile.common index e645d8457..d57d89614 100644 --- a/lammps/lib/colvars/Makefile.common +++ b/lammps/lib/colvars/Makefile.common @@ -39,6 +39,7 @@ COLVARS_SRCS = \ colvarcomp_distances.cpp \ colvarcomp_gpath.cpp \ colvarcomp_neuralnetwork.cpp \ + colvarcomp_torchann.cpp \ colvarcomp_combination.cpp \ colvarcomp_protein.cpp \ colvarcomp_rotations.cpp \ diff --git a/lammps/lib/colvars/Makefile.deps b/lammps/lib/colvars/Makefile.deps index dee9cbca6..c3029f45c 100644 --- a/lammps/lib/colvars/Makefile.deps +++ b/lammps/lib/colvars/Makefile.deps @@ -97,6 +97,11 @@ $(COLVARS_OBJ_DIR)colvarcomp_neuralnetwork.o: \ colvarproxy.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h \ colvar_neuralnetworkcompute.h +$(COLVARS_OBJ_DIR)colvarcomp_torchann.o: \ + colvarcomp_torchann.cpp colvarmodule.h colvars_version.h \ + colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \ + colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_io.h \ + colvarproxy_system.h colvarproxy_tcl.h $(COLVARS_OBJ_DIR)colvarcomp_combination.o: colvarcomp_combination.cpp \ colvarcomp.h colvarmodule.h colvars_version.h colvaratoms.h \ colvarproxy.h colvartypes.h ../../src/math_eigen_impl.h colvarvalue.h \ diff --git a/namd/colvars/Make.depends b/namd/colvars/Make.depends index af0333e7b..5598fd149 100644 --- a/namd/colvars/Make.depends +++ b/namd/colvars/Make.depends @@ -413,8 +413,27 @@ obj/colvarcomp_combination.o: \ colvars/src/colvardeps.h \ colvars/src/colvar.h \ colvars/src/colvar_arithmeticpath.h \ - colvars/src/colvar_geometricpath.h + colvars/src/colvar_geometricpath.h $(CXX) $(COLVARSCXXFLAGS) $(COPTO)obj/colvarcomp_combination.o $(COPTC) colvars/src/colvarcomp_combination.cpp +obj/colvarcomp_torchann.o: \ + obj/.exists \ + colvars/src/colvarcomp_torchann.cpp \ + colvars/src/colvarmodule.h \ + colvars/src/colvars_version.h \ + colvars/src/colvarvalue.h \ + colvars/src/colvartypes.h \ + colvars/src/colvarparse.h \ + colvars/src/colvarparams.h \ + colvars/src/colvar.h \ + colvars/src/colvardeps.h \ + colvars/src/colvarcomp.h \ + colvars/src/colvaratoms.h \ + colvars/src/colvarproxy.h \ + colvars/src/colvarproxy_tcl.h \ + colvars/src/colvarproxy_volmaps.h \ + colvars/src/colvar_arithmeticpath.h \ + colvars/src/colvar_geometricpath.h + $(CXX) $(COLVARSCXXFLAGS) $(COPTO)obj/colvarcomp_torchann.o $(COPTC) colvars/src/colvarcomp_torchann.cpp obj/colvarcomp_neuralnetwork.o: \ obj/.exists \ colvars/src/colvarcomp_neuralnetwork.cpp \ diff --git a/namd/colvars/src/Makefile.namd b/namd/colvars/src/Makefile.namd index 636f29681..c85f36f3b 100644 --- a/namd/colvars/src/Makefile.namd +++ b/namd/colvars/src/Makefile.namd @@ -20,6 +20,7 @@ COLVARSLIB = \ $(DSTDIR)/colvarcomp_volmaps.o \ $(DSTDIR)/colvarcomp_combination.o \ $(DSTDIR)/colvarcomp_neuralnetwork.o \ + $(DSTDIR)/colvarcomp_torchann.o \ $(DSTDIR)/colvar_neuralnetworkcompute.o \ $(DSTDIR)/colvardeps.o \ $(DSTDIR)/colvargrid.o \ diff --git a/namd/config.patch b/namd/config.patch new file mode 100644 index 000000000..90ef44065 --- /dev/null +++ b/namd/config.patch @@ -0,0 +1,59 @@ +diff --git a/config b/config +index 2bbb3e9a..d53d6d53 100755 +--- a/config ++++ b/config +@@ -56,6 +56,8 @@ function error_syntax { + echo ' --cuda-dlink arch=,code= (for cuFFT, may be repeated)' + echo ' --with-cuda-profiling (enables CUDA profiling with NVTX)' + echo ' --with-rocm-profiling (enables ROCm profiling with ROCtracer API)' ++ echo ' --with-colvars-torch (enables torchann in Colvars)' ++ echo ' --torch-prefix ' + echo '' + if [ -n "${PRINT_ARCH_LIST+set}" ]; then + ARCH_PAT='' +@@ -151,6 +153,7 @@ function error_exists { + use_mkl=0 + use_cuda=0 + use_hip=0 ++ use_colvars_torch=0 + use_cuda_prof=0 + use_rocm_prof=0 + use_memopt=0 +@@ -254,6 +257,19 @@ function error_exists { + ARCH_SUFFIX_ARG=$ARCH_SUFFIX_ARG-$1 + ;; + ++ --with-colvars-torch) ++ use_colvars_torch=1 ++ ;; ++ ++ --torch-prefix) ++ shift ++ if [ ! -d "$1" ]; then ++ echo "ERROR: No such directory $1" ++ error_syntax ++ fi ++ TORCH_PREFIX=$1 ++ ;; ++ + --with-debug) + use_debug=1 + ;; +@@ -1007,7 +1023,16 @@ function error_exists { + if [ -n "$CC_OPTS" ]; then + echo "COPTS = $CC_OPTS" >> Make.config + fi +- ++ if (( "$use_colvars_torch" )); then ++ echo "TORCHDIR = $TORCH_PREFIX" >> Make.config ++ echo 'TORCHINCFLAGS = $(COPTI)$(TORCHDIR)/include $(COPTI)$(TORCHDIR)/include/torch/csrc/api/include' >> Make.config ++ echo 'EXTRACOLVARSFLAGS = -std=c++17 -DTORCH $(TORCHINCFLAGS)' >> Make.config ++ if [[ $use_cuda && -f "${TORCH_PREFIX}/lib/libtorch_cuda.so" ]]; then ++ echo 'EXTRALINKLIBS = -Wl,-rpath,$(TORCHDIR)/lib -L$(TORCHDIR)/lib -ltorch -ltorch_cpu -lc10 -ltorch_cuda -lc10_cuda' >> Make.config ++ else ++ echo 'EXTRALINKLIBS = -Wl,-rpath,$(TORCHDIR)/lib -L$(TORCHDIR)/lib -ltorch -ltorch_cpu -lc10' >> Make.config ++ fi ++ fi + if (( $use_debug )); then + echo 'CXXOPTS = -g' >> Make.config + echo 'CXXTHREADOPTS = -g' >> Make.config diff --git a/src/Makefile b/src/Makefile index 96ba79d22..eebc27403 100644 --- a/src/Makefile +++ b/src/Makefile @@ -62,7 +62,7 @@ COLVARS_OBJS := $(COLVARS_SRCS:.cpp=.o) $(LEPTON_OBJS) $(CXX) $(CXXFLAGS) $(COLVARS_INCFLAGS) $(LEPTON_INCFLAGS) -c -o $@ $< $(COLVARS_LIB): Makefile.deps $(COLVARS_OBJS) - $(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS) $(LEPTON_OBJS) + $(AR) $(ARFLAGS) $(COLVARS_LIB) $(COLVARS_OBJS) $(LEPTON_OBJS) Makefile.deps: $(COLVARS_SRCS) @echo > $@ diff --git a/src/colvar.cpp b/src/colvar.cpp index 610ff9c85..d90089ce9 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -907,7 +907,9 @@ int colvar::init_components(std::string const &conf) error_code |= init_components_type(conf, "CV with support of the lepton custom function", "customColvar"); #endif error_code |= init_components_type(conf, "neural network CV for other CVs", "NeuralNetwork"); - +#ifdef TORCH + error_code |= init_components_type(conf, "CV defined by PyTorch artifical neural network models", "torchANN"); +#endif error_code |= init_components_type(conf, "total value of atomic map", "mapTotal"); // iterate over all available CVC in the map diff --git a/src/colvar.h b/src/colvar.h index 4dee696ce..00231984d 100644 --- a/src/colvar.h +++ b/src/colvar.h @@ -632,6 +632,7 @@ class colvar : public colvarparse, public colvardeps { class euler_psi; class euler_theta; class neuralNetwork; + class torchANN; class customColvar; // non-scalar components diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 1d7f2e839..19b24aeb4 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -28,6 +28,11 @@ #include "colvar.h" #include "colvar_geometricpath.h" +#ifdef TORCH +#include +#include +#endif + /// \brief Colvar component (base class for collective variables) /// @@ -1796,6 +1801,50 @@ class colvar::neuralNetwork virtual void apply_force(colvarvalue const &force); }; +#ifdef TORCH +// only when LibTorch is available +class colvar::torchANN + : public colvar::linearCombination +{ +protected: + torch::jit::script::Module nn; + /// the index of nn output component + size_t m_output_index; + bool use_double_input; + bool use_gpu; + // 1d tensor, concatenation of values of sub-cvcs + torch::Tensor input_tensor; + torch::Tensor nn_outputs; + torch::Tensor input_grad; + // record the initial index of of sub-cvcs in input_tensor + std::vector cvc_indices; +public: + torchANN(std::string const &conf); + virtual ~torchANN(); + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + + /// Redefined to handle periodicity + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to handle periodicity + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to handle periodicity + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to handle periodicity + virtual void wrap(colvarvalue &x_unwrapped) const; +}; +#else +class colvar::torchANN + : public colvar::componentDisabled +{ +public: + torchANN(std::string const &conf) : componentDisabled(conf) {} +}; +#endif // TORCH checking // \brief Colvar component: total value of a scalar map // (usually implemented as a grid by the simulation engine) diff --git a/src/colvarcomp_torchann.cpp b/src/colvarcomp_torchann.cpp new file mode 100644 index 000000000..2e7c03d0b --- /dev/null +++ b/src/colvarcomp_torchann.cpp @@ -0,0 +1,231 @@ +// -*- c++ -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. + +#ifdef TORCH + +#include "colvarmodule.h" +#include "colvar.h" +#include "colvarvalue.h" +#include "colvarparse.h" +#include "colvarcomp.h" + +colvar::torchANN::torchANN(std::string const &conf): linearCombination(conf) { + set_function_type("torchANN"); + + x.type(colvarvalue::type_scalar); + enable(f_cvc_scalar); + + if (period != 0.0) { + enable(f_cvc_periodic); + } + + if ((wrap_center != 0.0) && !is_enabled(f_cvc_periodic)) { + cvm::error("Error: wrapAround was defined in a torchANN component," + " but its period has not been set.\n"); + return; + } + + std::string model_file ; + get_keyval(conf, "modelFile", model_file, std::string("")); + try { + nn = torch::jit::load(model_file); + cvm::log("torch model loaded.") ; + } catch (const std::exception & e) { + cvm::error("Error: couldn't load libtorch model (see below).\n" + cvm::to_str(e.what())); + } + get_keyval(conf, "m_output_index", m_output_index, 0); + get_keyval(conf, "doubleInputTensor", use_double_input, false); + get_keyval(conf, "useGPU", use_gpu, false); + + cvc_indices.resize(cv.size(),0); + + size_t num_inputs = 0; + // compute total number of inputs of neural network + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) + { + num_inputs += cv[i_cv]->value().size() ; + if (i_cv < cv.size() - 1) + cvc_indices[i_cv+1] = num_inputs; + } + cvm::log("Input dimension of model: " + cvm::to_str(num_inputs)); + + // initialize the input tensor + auto options = torch::TensorOptions().dtype(torch::kFloat32).requires_grad(true); + + if (use_gpu) { + if (torch::cuda::is_available()) { + try { + nn.to(torch::kCUDA); + } catch(const std::exception & e) { + cvm::error("Failed to move model to GPU."); + use_gpu = false; + } + } else { + use_gpu = false; + cvm::log("GPU not available."); + } + } + + if (use_gpu) { + options = options.device(torch::kCUDA); + cvm::log("Use GPU."); + if (use_double_input) { + cvm::log("Data type reset to Float for GPU computation!"); + use_double_input = false; + } + } else { + cvm::log("Use CPU."); + } + + if (use_double_input) { // set type to double + options = options.dtype(torch::kFloat64); + nn.to(torch::kFloat64); + cvm::log("Model's dtype: kFloat64."); + } else { + cvm::log("Model's dtype: kFloat32."); + } + + input_tensor = torch::zeros({1,(long int) num_inputs}, options); + + try { // test the model + std::vector inputs={input_tensor}; + nn_outputs = nn.forward(inputs).toTensor()[0][m_output_index]; + cvm::log("Evaluating model with zero tensor succeeded."); + } catch (const std::exception & e) { + cvm::error("Error: evaluating libtorch model with zero tensor failed (see below).\n" + cvm::to_str(e.what())); + } +} + +colvar::torchANN::~torchANN() { +} + +void colvar::torchANN::calc_value() { + + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) + cv[i_cv]->calc_value(); + + if (use_gpu) + input_tensor = input_tensor.to(torch::kCPU); + + // set input tensor with no_grad + { + torch::NoGradGuard no_grad; + size_t l = 0; + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + const colvarvalue& current_cv_value = cv[i_cv]->value(); + if (current_cv_value.type() == colvarvalue::type_scalar) { + input_tensor[0][l++] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)); + } else { + for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) + input_tensor[0][l++] = cv[i_cv]->sup_coeff * current_cv_value[j_elem]; + } + } + } + + if (use_gpu) + input_tensor = input_tensor.to(torch::kCUDA); + + std::vector inputs={input_tensor}; + + // evaluate the value of function + nn_outputs = nn.forward(inputs).toTensor()[0][m_output_index]; + + input_grad = torch::autograd::grad({nn_outputs}, {input_tensor})[0][0]; + + if (use_gpu) + input_grad = input_grad.to(torch::kCPU); + + x = nn_outputs.item() ; + + this->wrap(x); +} + +void colvar::torchANN::calc_gradients() { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_gradients(); + if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + // get the initial index of this cvc + size_t l = cvc_indices[i_cv]; + for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) { + // get derivative of neural network wrt its input + const cvm::real factor = input_grad[l+j_elem].item(); + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * factor * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; + } + } + } + } + } +} + +void colvar::torchANN::apply_force(colvarvalue const &force) { + + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + // If this CV us explicit gradients, then atomic gradients is already calculated + // We can apply the force to atom groups directly + if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); + } + } else { + const colvarvalue& current_cv_value = cv[i_cv]->value(); + colvarvalue cv_force(current_cv_value); + cv_force.reset(); + const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + // get the initial index of this cvc + size_t l = cvc_indices[i_cv]; + for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { + cv_force[j_elem] = factor_polynomial * input_grad[l+j_elem].item() * force.real_value; + } + cv[i_cv]->apply_force(cv_force); + } + } +} + +cvm::real colvar::torchANN::dist2(colvarvalue const &x1, colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + if (is_enabled(f_cvc_periodic)) + diff = (diff < - period * 0.5 ? diff + period : (diff > period * 0.5 ? diff - period : diff)); + return diff * diff; +} + +colvarvalue colvar::torchANN::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + if (is_enabled(f_cvc_periodic)) + diff = (diff < - period * 0.5 ? diff + period : (diff > period * 0.5 ? diff - period : diff)); + return 2.0 * diff; +} + + +colvarvalue colvar::torchANN::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + if (is_enabled(f_cvc_periodic)) + diff = (diff < - period * 0.5 ? diff + period : (diff > period * 0.5 ? diff - period : diff)); + return (-2.0) * diff; +} + + +void colvar::torchANN::wrap(colvarvalue &x_unwrapped) const +{ + if (!is_enabled(f_cvc_periodic)) { + return; + } + cvm::real shift = + cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5); + x_unwrapped.real_value -= shift * period; +} + +#endif diff --git a/src/colvarmodule_refs.h b/src/colvarmodule_refs.h index b3f584242..74aad5144 100644 --- a/src/colvarmodule_refs.h +++ b/src/colvarmodule_refs.h @@ -587,3 +587,6 @@ feature_count_[std::string("Scripted functions (Tcl)")] = 0; feature_paper_map_[std::string("Scripted functions (Tcl)")] = "n/a"; + + feature_count_[std::string("torchANN colvar component")] = 0; + feature_paper_map_[std::string("torchANN colvar component")] = "n/a"; diff --git a/update-colvars-code.sh b/update-colvars-code.sh index 3ddcb4880..c75ab3083 100755 --- a/update-colvars-code.sh +++ b/update-colvars-code.sh @@ -317,6 +317,11 @@ then condcopy "${src}" "${target}/lib/colvars/${tgt}" done + # Update cmake file + src=${source}/lammps/COLVARS.cmake + tgt=$(basename ${src}) + condcopy "${source}/lammps/COLVARS.cmake" "${target}/cmake/Modules/Packages/${tgt}" + for src in \ ${source}/lammps/src/COLVARS/colvarproxy_lammps{.cpp,.h,_version.h} \ ${source}/lammps/src/COLVARS/fix_colvars.{cpp,h} \ @@ -384,6 +389,10 @@ then patch -p1 -N -d ${target} < namd/Makefile.patch fi + # For torchann + patch -F 12 -p1 -N -d ${target} < namd/config.patch + sed -i '/^COLVARSCXXFLAGS =*/s/$/ $(EXTRACOLVARSFLAGS)/' ${target}/Makefile + # Copy library files to the "colvars" folder for src in ${source}/src/*.h ${source}/src/*.cpp do \ diff --git a/vmd/src/colvars_files.pl b/vmd/src/colvars_files.pl index b6bfdbf2a..a79df23f1 100644 --- a/vmd/src/colvars_files.pl +++ b/vmd/src/colvars_files.pl @@ -22,6 +22,7 @@ 'colvar.C', 'colvarcomp_neuralnetwork.C', 'colvar_neuralnetworkcompute.C', + 'colvarcomp_torchann.C', 'colvarcomp.C', 'colvarcomp_alchlambda.C', 'colvarcomp_angles.C',