From ddb743b124b660014ae7a1cb2749459108eabc2d Mon Sep 17 00:00:00 2001 From: Hubert Date: Fri, 28 Jul 2023 12:45:17 +0200 Subject: [PATCH 1/5] support GROMACS 2023 --- .github/workflows/test-backends.yml | 12 + gromacs/gromacs-2023.x.patch | 514 ++++++++++++++++ gromacs/gromacs-2023.x/CMakeLists.txt | 35 ++ .../gromacs-2023.x/colvarproxy_gromacs.cpp | 575 ++++++++++++++++++ gromacs/gromacs-2023.x/colvarproxy_gromacs.h | 144 +++++ update-colvars-code.sh | 17 +- 6 files changed, 1293 insertions(+), 4 deletions(-) create mode 100644 gromacs/gromacs-2023.x.patch create mode 100644 gromacs/gromacs-2023.x/CMakeLists.txt create mode 100644 gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp create mode 100644 gromacs/gromacs-2023.x/colvarproxy_gromacs.h diff --git a/.github/workflows/test-backends.yml b/.github/workflows/test-backends.yml index 266972893..edaf7caff 100644 --- a/.github/workflows/test-backends.yml +++ b/.github/workflows/test-backends.yml @@ -138,3 +138,15 @@ jobs: path_compile_script: devel-tools/compile-gromacs.sh test_lib_directory: gromacs/tests/library rpath_exe: install/bin/gmx_d + + gromacs-2023: + name: GROMACS 2023 + if: github.event_name == 'pull_request' || contains(github.event.head_commit.message, 'test-gromacs-2023') + uses: ./.github/workflows/backend-template.yml + with: + backend_name: GROMACS-2023 + backend_repo: gromacs/gromacs + backend_repo_ref: release-2023 + path_compile_script: devel-tools/compile-gromacs.sh + test_lib_directory: gromacs/tests/library + rpath_exe: install/bin/gmx_d \ No newline at end of file diff --git a/gromacs/gromacs-2023.x.patch b/gromacs/gromacs-2023.x.patch new file mode 100644 index 000000000..ef5a373b6 --- /dev/null +++ b/gromacs/gromacs-2023.x.patch @@ -0,0 +1,514 @@ +diff --git a/CMakeLists.txt b/CMakeLists.txt +index 60636ec..0cfab4b 100644 +--- a/CMakeLists.txt ++++ b/CMakeLists.txt +@@ -634,6 +634,10 @@ include(gmxManageLmfit) + + include(gmxManageMuparser) + ++include(gmxManageColvars) ++ ++include(gmxManageLepton) ++ + ################################################## + # Process SIMD instruction settings + ################################################## +diff --git a/api/legacy/include/gromacs/mdtypes/inputrec.h b/api/legacy/include/gromacs/mdtypes/inputrec.h +index 8ddc9c1..21b3221 100644 +--- a/api/legacy/include/gromacs/mdtypes/inputrec.h ++++ b/api/legacy/include/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; +@@ -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; ++ ++ /* COLVARS */ ++ bool bColvars = false; /* Do we do colvars calculations ? */ ++ colvarproxy_gromacs *colvars_proxy = nullptr; /* The object for the colvars calculations */ + }; + + int tcouple_min_integration_steps(TemperatureCoupling etc); +diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt +index 114bfd9..bb03eb4 100644 +--- a/src/gromacs/CMakeLists.txt ++++ b/src/gromacs/CMakeLists.txt +@@ -146,6 +146,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) ++list(APPEND libgromacs_object_library_dependencies lepton) ++ + # 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() + add_library(libgromacs ${LIBGROMACS_SOURCES}) + endif() + ++gmx_include_colvars_headers() ++ + 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 3c4987f..53e6b91 100644 +--- a/src/gromacs/applied_forces/CMakeLists.txt ++++ b/src/gromacs/applied_forces/CMakeLists.txt +@@ -67,6 +67,7 @@ gmx_add_libgromacs_sources( + add_subdirectory(awh) + 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 2f2cb69..59611d2 100644 +--- a/src/gromacs/fileio/checkpoint.cpp ++++ b/src/gromacs/fileio/checkpoint.cpp +@@ -69,6 +69,7 @@ + #include "gromacs/mdtypes/pullhistory.h" + #include "gromacs/mdtypes/state.h" + #include "gromacs/mdtypes/swaphistory.h" ++#include "gromacs/mdtypes/colvarshistory.h" + #include "gromacs/modularsimulator/modularsimulator.h" + #include "gromacs/trajectory/trajectoryframe.h" + #include "gromacs/utility/arrayref.h" +@@ -1289,6 +1290,14 @@ static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderC + { + contents->isModularSimulatorCheckpoint = false; + } ++ if (contents->file_version >= CheckPointVersion::Colvars) ++ { ++ do_cpt_int_err(xd, "colvars atoms", &contents->ecolvars, list); ++ } ++ else ++ { ++ contents->ecolvars = 0; ++ } + } + + static int do_cpt_footer(XDR* xd, CheckPointVersion file_version) +@@ -2041,6 +2050,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) ++{ ++ ++ if (ecolvars == 0) ++ { ++ return 0; ++ } ++ ++ colvarsstate->bFromCpt = bRead; ++ colvarsstate->n_atoms = ecolvars; ++ ++ /* Write data */ ++ char buf[STRLEN]; ++ sprintf(buf, "Colvars atoms in reference structure : %d", ecolvars); ++ sprintf(buf, "Colvars xa_old"); ++ if (bRead) ++ { ++ snew(colvarsstate->xa_old_whole, colvarsstate->n_atoms); ++ do_cpt_n_rvecs_err(xd, buf, colvarsstate->n_atoms, colvarsstate->xa_old_whole, list); ++ } ++ else ++ { ++ do_cpt_n_rvecs_err(xd, buf, colvarsstate->n_atoms, colvarsstate->xa_old_whole_p, list); ++ } ++ ++ return 0; ++} ++ ++ + static int do_cpt_correlation_grid(XDR* xd, + gmx_bool bRead, + gmx_unused int fflags, +@@ -2453,6 +2492,7 @@ void write_checkpoint_data(t_fileio* fp, + observablesHistory->swapHistory.get(), + nullptr) + < 0) ++ || (do_cpt_colvars(gmx_fio_getxdr(fp), FALSE, headerContents.ecolvars, observablesHistory->colvarsHistory.get(), nullptr) < 0) + || (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?"); +@@ -2853,6 +2893,18 @@ static void read_checkpoint(const std::filesystem::path& fn, + cp_error(); + } + ++ if (headerContents->ecolvars != 0 && observablesHistory->colvarsHistory == nullptr) ++ { ++ observablesHistory->colvarsHistory = std::make_unique(colvarshistory_t{}); ++ } ++ ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, headerContents->ecolvars, ++ observablesHistory->colvarsHistory.get(), nullptr); ++ if (ret) ++ { ++ cp_error(); ++ } ++ ++ + std::vector outputfiles; + ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version); + if (ret) +@@ -3032,6 +3084,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) ++ { ++ cp_error(); ++ } ++ + ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version); + + if (ret) +@@ -3157,6 +3216,12 @@ void list_checkpoint(const std::filesystem::path& fn, FILE* out) + ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out); + } + ++ if (ret == 0) ++ { ++ colvarshistory_t colvarshist = {}; ++ ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, headerContents.ecolvars, &colvarshist, out); ++ } ++ + if (ret == 0) + { + std::vector outputfiles; +diff --git a/src/gromacs/fileio/checkpoint.h b/src/gromacs/fileio/checkpoint.h +index 4811cb8..d2a32d0 100644 +--- a/src/gromacs/fileio/checkpoint.h ++++ b/src/gromacs/fileio/checkpoint.h +@@ -220,6 +220,8 @@ enum class CheckPointVersion : int + ModularSimulator, + //! Added local (per walker) weight contribution to each point in AWH. + AwhLocalWeightSum, ++ //! Added Colvars ++ Colvars, + //! The total number of checkpoint versions. + Count, + //! Current version +@@ -302,6 +304,8 @@ struct CheckpointHeaderContents + int nED; + //! Enum for coordinate swapping. + SwapType eSwapCoords; ++ //! Colvars ++ int ecolvars; + //! Whether the checkpoint was written by modular simulator. + bool isModularSimulatorCheckpoint = false; + }; +diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp +index 6fbd423..c57b4db 100644 +--- a/src/gromacs/mdlib/energyoutput.cpp ++++ b/src/gromacs/mdlib/energyoutput.cpp +@@ -255,7 +255,7 @@ EnergyOutput::EnergyOutput(ener_file* fp_ene, + bEner_[F_DISPCORR] = (inputrec.eDispCorr != DispersionCorrectionType::No); + bEner_[F_DISRESVIOL] = (gmx_mtop_ftype_count(mtop, F_DISRES) > 0); + 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 +index 3ea5b73..c88bdc6 100644 +--- a/src/gromacs/mdlib/mdoutf.cpp ++++ b/src/gromacs/mdlib/mdoutf.cpp +@@ -60,6 +60,7 @@ + #include "gromacs/mdtypes/observableshistory.h" + #include "gromacs/mdtypes/state.h" + #include "gromacs/mdtypes/swaphistory.h" ++#include "gromacs/mdtypes/colvarshistory.h" + #include "gromacs/timing/wallcycle.h" + #include "gromacs/topology/topology.h" + #include "gromacs/utility/baseversion.h" +@@ -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); ++ + CheckpointHeaderContents headerContents = { CheckPointVersion::UnknownVersion0, + { 0 }, + { 0 }, +diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp +index a839ab7..625b437 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; +@@ -647,6 +649,16 @@ static void computeSpecialForces(FILE* fplog, + */ + if (stepWork.computeForces) + { ++ ++ /* COLVARS */ ++ /* Colvars Module needs some updated data - just PBC & step number - before calling its ForceProvider */ ++ if (inputrec.bColvars) ++ { ++ t_pbc pbc; ++ set_pbc(&pbc, inputrec.pbcType, box); ++ inputrec.colvars_proxy->update_data(cr, step, pbc, box, didNeighborSearch); ++ } ++ + gmx::ForceProviderInput forceProviderInput( + x, + mdatoms->homenr, +diff --git a/src/gromacs/mdrun/legacymdrunoptions.h b/src/gromacs/mdrun/legacymdrunoptions.h +index 20d793e..2d9b27b 100644 +--- a/src/gromacs/mdrun/legacymdrunoptions.h ++++ b/src/gromacs/mdrun/legacymdrunoptions.h +@@ -124,7 +124,10 @@ public: + { efTOP, "-mp", "membed", ffOPTRD }, + { efNDX, "-mn", "membed", ffOPTRD }, + { efXVG, "-if", "imdforces", ffOPTWR }, +- { efXVG, "-swap", "swapions", ffOPTWR } } }; ++ { 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 +index 6da922e..89088f3 100644 +--- a/src/gromacs/mdrun/replicaexchange.cpp ++++ b/src/gromacs/mdrun/replicaexchange.cpp +@@ -624,6 +624,7 @@ static void exchange_state(const gmx_multisim_t* ms, int b, t_state* state) + exchange_doubles(ms, b, &state->baros_integral, 1); + exchange_rvecs(ms, b, state->x.rvec_array(), state->natoms); + 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 acadeab..faf31b5 100644 +--- a/src/gromacs/mdrun/runner.cpp ++++ b/src/gromacs/mdrun/runner.cpp +@@ -131,6 +131,7 @@ + #include "gromacs/mdtypes/mdrunoptions.h" + #include "gromacs/mdtypes/multipletimestepping.h" + #include "gromacs/mdtypes/observableshistory.h" ++#include "gromacs/mdtypes/colvarshistory.h" + #include "gromacs/mdtypes/observablesreducer.h" + #include "gromacs/mdtypes/simulation_workload.h" + #include "gromacs/mdtypes/state.h" +@@ -175,6 +176,8 @@ + #include "gromacs/utility/smalloc.h" + #include "gromacs/utility/stringutil.h" + ++#include "gromacs/applied_forces/colvars/colvarproxy_gromacs.h" ++ + #include "isimulator.h" + #include "membedholder.h" + #include "replicaexchange.h" +@@ -2142,6 +2145,51 @@ int Mdrunner::mdrunner() + mdrunOptions.imdOptions, + startingBehavior); + ++ /* COLVARS */ ++ if (opt2bSet("-colvars",filenames.size(), filenames.data())) ++ { ++ ++ gmx::ArrayRef filenames_colvars; ++ std::string filename_restart; ++ std::string prefix; ++ ++ inputrec->bColvars = TRUE; ++ ++ /* Retrieve filenames */ ++ filenames_colvars = opt2fns("-colvars", filenames.size(), filenames.data()); ++ if (opt2bSet("-colvars_restart",filenames.size(), filenames.data())) ++ { ++ filename_restart = opt2fn("-colvars_restart",filenames.size(), filenames.data()); ++ } ++ ++ /* Determine the prefix for the colvars output files, based on the logfile name. */ ++ std::string logfile = ftp2fn(efLOG, filenames.size(), filenames.data()); ++ /* 4 = ".log".length() */ ++ if(logfile.length() > 4) ++ { ++ prefix = logfile.substr(0,logfile.length()-4); ++ } ++ ++ inputrec->colvars_proxy = new colvarproxy_gromacs(); ++ inputrec->colvars_proxy->init(inputrec.get(),inputrec->init_step,mtop, &observablesHistory, prefix, filenames_colvars,filename_restart, cr, MAIN(cr) ? globalState->x.rvec_array() : nullptr, ++ MAIN(cr) ? &globalState->xa_old_whole_colvars : nullptr, MAIN(cr) ? &globalState->n_colvars_atoms : nullptr); ++ fr->forceProviders->addForceProvider(inputrec->colvars_proxy); ++ } ++ else ++ { ++ inputrec->bColvars = FALSE; ++ if (opt2bSet("-colvars_restart",filenames.size(), filenames.data())) ++ { ++ gmx_fatal(FARGS, "-colvars_restart can only be used together with the -colvars option."); ++ } ++ if(observablesHistory.colvarsHistory) { ++ gmx_fatal(FARGS, ++ "The checkpoint is from a run with colvars, " ++ "but the current run did not specify the -colvars option. " ++ "Either specify the -colvars option to mdrun, or do not use this checkpoint file."); ++ } ++ } ++ + if (haveDDAtomOrdering(*cr)) + { + GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP"); +@@ -2336,6 +2384,16 @@ int Mdrunner::mdrunner() + releaseDevice(deviceInfo); + } + ++ /* COLVARS */ ++ if (inputrec->bColvars) ++ { ++ GMX_RELEASE_ASSERT(inputrec->colvars_proxy, "inputrec->colvars_proxy was NULL while colvars module was enabled."); ++ ++ inputrec->colvars_proxy->finish(cr); ++ delete inputrec->colvars_proxy; ++ inputrec->colvars_proxy = nullptr; ++ } ++ + /* 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 +new file mode 100644 +index 0000000..6605e6f +--- /dev/null ++++ b/src/gromacs/mdtypes/colvarshistory.h +@@ -0,0 +1,20 @@ ++#ifndef GMX_MDLIB_COLVARSHISTORY_H ++#define GMX_MDLIB_COLVARSHISTORY_H ++ ++#include "gromacs/math/vectypes.h" ++#include "gromacs/utility/basedefinitions.h" ++ ++/* Helper structure to be able to make colvars group(s) whole ++ * ++ * To also make colvars group(s) whole, we save the last whole configuration ++ * of the atoms in the checkpoint file. ++ */ ++typedef struct colvarshistory_t ++{ ++ gmx_bool bFromCpt; // Did we start from a checkpoint file? ++ int n_atoms; // Number of colvars atoms ++ rvec* xa_old_whole; // Last known whole positions of the colvars atoms ++ rvec* xa_old_whole_p; // Pointer to these positions ++} colvarshistory_t; ++ ++#endif +diff --git a/src/gromacs/mdtypes/observableshistory.cpp b/src/gromacs/mdtypes/observableshistory.cpp +index 12230ee..13f5df1 100644 +--- a/src/gromacs/mdtypes/observableshistory.cpp ++++ b/src/gromacs/mdtypes/observableshistory.cpp +@@ -40,6 +40,7 @@ + #include "gromacs/mdtypes/energyhistory.h" + #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 +index 6ed0e3f..172ead0 100644 +--- a/src/gromacs/mdtypes/observableshistory.h ++++ b/src/gromacs/mdtypes/observableshistory.h +@@ -58,6 +58,7 @@ class energyhistory_t; + class PullHistory; + 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 6a38311..872673f 100644 +--- a/src/gromacs/mdtypes/state.cpp ++++ b/src/gromacs/mdtypes/state.cpp +@@ -370,7 +370,9 @@ t_state::t_state() : + dfhist(nullptr), + awhHistory(nullptr), + ddp_count(0), +- ddp_count_cg_gl(0) ++ ddp_count_cg_gl(0), ++ xa_old_whole_colvars(nullptr), ++ n_colvars_atoms(0) + + { + clear_mat(box); +diff --git a/src/gromacs/mdtypes/state.h b/src/gromacs/mdtypes/state.h +index c8c05e8..2026f26 100644 +--- a/src/gromacs/mdtypes/state.h ++++ b/src/gromacs/mdtypes/state.h +@@ -285,6 +285,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 +--- a/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml ++++ b/src/programs/mdrun/tests/refdata/MdrunTest_WritesHelp.xml +@@ -6,7 +6,8 @@ + gmx [-s [<.tpr>]] [-cpi [<.cpt>]] [-table [<.xvg>]] [-tablep [<.xvg>]] + [-tableb [<.xvg> [...]]] [-rerun [<.xtc/.trr/...>]] [-ei [<.edi>]] + [-multidir [<dir> [...]]] [-awh [<.xvg>]] [-membed [<.dat>]] +- [-mp [<.top>]] [-mn [<.ndx>]] [-o [<.trr/.cpt/...>]] [-x [<.xtc/.tng>]] ++ [-mp [<.top>]] [-mn [<.ndx>]] [-colvars [<.dat> [...]]] ++ [-colvars_restart [<.dat>]] [-o [<.trr/.cpt/...>]] [-x [<.xtc/.tng>]] + [-cpo [<.cpt>]] [-c [<.gro/.g96/...>]] [-e [<.edr>]] [-g [<.log>]] + [-dhdl [<.xvg>]] [-field [<.xvg>]] [-tpi [<.xvg>]] [-tpid [<.xvg>]] + [-eo [<.xvg>]] [-px [<.xvg>]] [-pf [<.xvg>]] [-ro [<.xvg>]] +@@ -163,6 +164,10 @@ Options to specify input files: + Topology file + -mn [<.ndx>] (membed.ndx) (Opt.) + Index file ++ -colvars [<.dat> [...]] (colvars.dat) (Opt.) ++ Generic data file ++ -colvars_restart [<.dat>] (colvars.dat) (Opt.) ++ Generic data file + + Options to specify output files: + diff --git a/gromacs/gromacs-2023.x/CMakeLists.txt b/gromacs/gromacs-2023.x/CMakeLists.txt new file mode 100644 index 000000000..551c252d3 --- /dev/null +++ b/gromacs/gromacs-2023.x/CMakeLists.txt @@ -0,0 +1,35 @@ +# +# This file is part of the GROMACS molecular simulation package. +# +# Copyright 2014- The GROMACS Authors +# and the project initiators Erik Lindahl, Berk Hess and David van der Spoel. +# Consult the AUTHORS/COPYING files and https://www.gromacs.org for details. +# +# GROMACS is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public License +# as published by the Free Software Foundation; either version 2.1 +# of the License, or (at your option) any later version. +# +# GROMACS is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with GROMACS; if not, see +# https://www.gnu.org/licenses, or write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +# +# If you want to redistribute modifications to GROMACS, please +# consider that scientific software is very special. Version +# control is crucial - bugs must be traceable. We will be happy to +# consider code for inclusion in the official distribution, but +# derived work must not be called official GROMACS. Details are found +# in the README & COPYING files - if they are missing, get the +# official version at https://www.gromacs.org. +# +# To help us fund GROMACS development, we humbly ask that you cite +# the research papers on the package. Check out https://www.gromacs.org. + +gmx_add_libgromacs_sources(colvarproxy_gromacs.cpp) + diff --git a/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp b/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp new file mode 100644 index 000000000..59b8661f5 --- /dev/null +++ b/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp @@ -0,0 +1,575 @@ +/// -*- c++ -*- + +#include +#include +#include +#include + + +#include "gromacs/math/units.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/forceoutput.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/futil.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/domdec/domdec_struct.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/domdec/ga2la.h" +#include "gromacs/mdtypes/colvarshistory.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/mdlib/broadcaststructs.h" + + +#include "colvarproxy_gromacs.h" + + +//************************************************************ +// colvarproxy_gromacs +colvarproxy_gromacs::colvarproxy_gromacs() : colvarproxy() {} + +// Colvars Initialization +void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &mtop, + ObservablesHistory* oh, + const std::string &prefix, + gmx::ArrayRef filenames_config, + const std::string &filename_restart, + const t_commrec *cr, + const rvec x[], + rvec **xa_old_whole_colvars_state_p, + int *n_colvars_atoms_state_p) { + + + // Initialize colvars. + first_timestep = true; + restart_frequency_s = 0; + + // User-scripted forces are not available in GROMACS + have_scripts = false; + + angstrom_value_ = 0.1; + + boltzmann_ = gmx::c_boltz; + + // Get the thermostat temperature. + // NOTE: Considers only the first temperature coupling group! + set_target_temperature(ir->opts.ref_t[0]); + + // GROMACS random number generation. + // Seed with the mdp parameter ld_seed, the Langevin dynamics seed. + rng.seed(ir->ld_seed); + + /// Handle input filenames and prefix/suffix for colvars files. + /// + /// filename_config is the colvars configuration file collected from "-colvars" option. + /// The output prefix will be the prefix of Gromacs log filename. + /// or "output" otherwise. + /// + /// For restart, 'filename_restart' is the colvars input file for restart, + /// set by the "-cv_restart" option. It will be NULL otherwise. + /// + + if(!prefix.empty()) + { + output_prefix_str = prefix; + } + else { + output_prefix_str = "output"; + } + + restart_output_prefix_str = prefix + ".restart"; + + colvars_restart = false; + + if(!filename_restart.empty()) + { + colvars_restart = true; + input_prefix_str = filename_restart; + input_prefix_str.erase(input_prefix_str.rfind(".dat")); + input_prefix_str.erase(input_prefix_str.rfind(".colvars.state")); + } + + // Retrieve masses and charges from input file + updated_masses_ = updated_charges_ = true; + + // Get GROMACS timestep (picosecond to femtosecond) + set_integration_timestep(ir->delta_t * 1000.0); + // Retrieve the topology of all atoms + gmx_atoms = gmx_mtop_global_atoms(mtop); + + // Read configuration file and set up the proxy only on the master node. + if (MASTER(cr)) + { + + // initiate module: this object will be the communication proxy + // colvarmodule pointer is only defined on the Master due to the static pointer to colvarproxy. + colvars = new colvarmodule(this); + + version_int = get_version_from_string(COLVARPROXY_VERSION); + + colvars->cite_feature("GROMACS engine"); + colvars->cite_feature("Colvars-GROMACS interface"); + + if (cvm::debug()) { + log("Initializing the colvars proxy object.\n"); + } + + cvm::log("Using GROMACS interface, version "+ + cvm::to_str(COLVARPROXY_VERSION)+".\n"); + + auto i = filenames_config.begin(); + for(; i != filenames_config.end(); ++i) { + add_config("configfile", i->c_str()); + } + + colvarproxy::parse_module_config(); + colvars->update_engine_parameters(); + colvars->setup_input(); + + // Citation Reporter + cvm::log(std::string("\n")+colvars->feature_report(0)+std::string("\n")); + + colvars->setup_output(); + + if (step != 0) { + cvm::log("Initializing step number to "+cvm::to_str(step)+".\n"); + } + + colvars->it = colvars->it_restart = step; + + + } // end master + + + // MPI initialisation + + // Initialise attributs for the MPI communication + if(MASTER(cr)) { + // Retrieve the number of colvar atoms + n_colvars_atoms = atoms_ids.size(); + // Copy their global indices + ind = atoms_ids.data(); // This has to be updated if the vector is reallocated + } + + + if(PAR(cr)) { + // Let the other nodes know the number of colvar atoms. + block_bc(cr->mpi_comm_mygroup, n_colvars_atoms); + + // Initialise atoms_new_colvar_forces on non-master nodes + if(!MASTER(cr)) { + atoms_new_colvar_forces.reserve(n_colvars_atoms); + } + } + + snew(x_colvars_unwrapped, n_colvars_atoms); + snew(xa_ind, n_colvars_atoms); + snew(xa_shifts, n_colvars_atoms); + snew(xa_eshifts, n_colvars_atoms); + snew(xa_old_whole, n_colvars_atoms); + snew(f_colvars, n_colvars_atoms); + + // Prepare data + + // Manage restart with .cpt + if (MASTER(cr)) + { + /* colvarsHistory is the struct holding the data saved in the cpt + + If we dont start with from a .cpt, prepare the colvarsHistory struct for proper .cpt writing, + If we did start from .cpt, we copy over the last whole structures from .cpt, + In any case, for subsequent checkpoint writing, we set the pointers (xa_old_whole_p) in + the xa_old_whole arrays, which contain the correct PBC representation of + colvars atoms at the last time step. + */ + + if (oh->colvarsHistory == nullptr) + { + oh->colvarsHistory = std::make_unique(colvarshistory_t{}); + } + colvarshistory_t *colvarshist = oh->colvarsHistory.get(); + + + snew(colvarshist->xa_old_whole_p, n_colvars_atoms); + + /* We always need the last whole positions such that + * in the next time step we can make the colvars atoms whole again in PBC */ + if (colvarshist->bFromCpt) + { + for (int i = 0; i < n_colvars_atoms; i++) + { + copy_rvec(colvarshist->xa_old_whole[i], xa_old_whole[i]); + } + } + else + { + colvarshist->n_atoms = n_colvars_atoms; + for (int i = 0; i < n_colvars_atoms; i++) + { + int ii = ind[i]; + copy_rvec(x[ii], xa_old_whole[i]); + } + } + + /* For subsequent checkpoint writing, set the pointers (xa_old_whole_p) to the xa_old_whole + * arrays that get updated at every NS step */ + colvarshist->xa_old_whole_p = xa_old_whole; + //Initialize number of colvars atoms from the global state + *n_colvars_atoms_state_p = n_colvars_atoms; + // Point the shifts array from the global state to the local shifts array + *xa_old_whole_colvars_state_p = xa_old_whole; + } + + + // Communicate initial coordinates and global indices to all processes + if (PAR(cr)) + { + nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, xa_old_whole); + snew_bc(MASTER(cr), ind, n_colvars_atoms); + nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, ind); + } + + // Serial Run + if (!PAR(cr)) + { + nat_loc = n_colvars_atoms; + nalloc_loc = n_colvars_atoms; + ind_loc = ind; + + // xa_ind[i] needs to be set to i for serial runs + for (int i = 0; i < n_colvars_atoms; i++) + { + xa_ind[i] = i; + } + } + + if (MASTER(cr) && cvm::debug()) { + cvm::log ("atoms_ids = "+cvm::to_str (atoms_ids)+"\n"); + cvm::log ("atoms_refcount = "+cvm::to_str (atoms_refcount)+"\n"); + cvm::log ("positions = "+cvm::to_str (atoms_positions)+"\n"); + cvm::log ("atoms_new_colvar_forces = "+cvm::to_str (atoms_new_colvar_forces)+"\n"); + cvm::log (cvm::line_marker); + log("done initializing the colvars proxy object.\n"); + } + + +} // End colvars initialization. + + +colvarproxy_gromacs::~colvarproxy_gromacs() +{} + +void colvarproxy_gromacs::finish(const t_commrec *cr) +{ + if(MASTER(cr)) { + colvars->write_restart_file(output_prefix_str+".colvars.state"); + colvars->write_output_files(); + } +} + +cvm::real colvarproxy_gromacs::rand_gaussian() +{ + return normal_distribution(rng); +} + +size_t colvarproxy_gromacs::restart_frequency() +{ + return restart_frequency_s; +} + +// **************** PERIODIC BOUNDARY CONDITIONS **************** +// Get the PBC-aware distance vector between two positions +cvm::rvector colvarproxy_gromacs::position_distance (cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const +{ + rvec r1, r2, dr; + r1[0] = pos1.x; + r1[1] = pos1.y; + r1[2] = pos1.z; + r2[0] = pos2.x; + r2[1] = pos2.y; + r2[2] = pos2.z; + + pbc_dx(&gmx_pbc, r2, r1, dr); + return cvm::atom_pos( dr[0], dr[1], dr[2] ); +} + + +void colvarproxy_gromacs::log (std::string const &message) +{ + std::istringstream is(message); + std::string line; + while (std::getline(is, line)) + // Gromacs prints messages on the stderr FILE + fprintf(stderr, "colvars: %s\n", line.c_str()); +} + +void colvarproxy_gromacs::error (std::string const &message) +{ + // In GROMACS, all errors are fatal. + fatal_error (message); +} + +void colvarproxy_gromacs::fatal_error (std::string const &message) +{ + log(message); + if (!cvm::debug()) + log("If this error message is unclear, " + "try recompiling with -DCOLVARS_DEBUG.\n"); + gmx_fatal(FARGS,"Error in collective variables module.\n"); +} + +void colvarproxy_gromacs::exit (std::string const gmx_unused &message) +{ + gmx_fatal(FARGS,"SUCCESS: %s\n", message.c_str()); +} + +int colvarproxy_gromacs::load_atoms (char const gmx_unused *filename, std::vector gmx_unused &atoms, + std::string const gmx_unused &pdb_field, double const gmx_unused pdb_field_value) +{ + cvm::error("Selecting collective variable atoms " + "from a PDB file is currently not supported.\n"); + return COLVARS_NOT_IMPLEMENTED; +} + +int colvarproxy_gromacs::load_coords (char const gmx_unused *filename, std::vector gmx_unused &pos, + const std::vector gmx_unused &indices, std::string const gmx_unused &pdb_field_str, + double const gmx_unused pdb_field_value) +{ + cvm::error("Loading atoms coordinates from a PDB or GRO file is currently not supported." + "Please use an XYZ file.\n"); + return COLVARS_NOT_IMPLEMENTED; +} + +int colvarproxy_gromacs::set_unit_system(std::string const &units_in, bool /*colvars_defined*/) +{ + if (units_in != "gromacs") { + cvm::error("Specified unit system \"" + units_in + "\" is unsupported in Gromacs. Supported units are \"gromacs\" (nm, kJ/mol).\n"); + return COLVARS_ERROR; + } + return COLVARS_OK; +} + +int colvarproxy_gromacs::backup_file (char const *filename) +{ + // Incremental gromacs backup system will be use only for those file + if (std::string(filename).rfind(std::string(".colvars.traj")) != std::string::npos) { + + // GROMACS function + make_backup(filename); + + // Otherwise, just keep one backup. + } else { + + //Handle filename of the backup file + const char *extension = ".old"; + char *backup = new char[strlen(filename)+strlen(extension)+1]; + strcpy(backup, filename); + strcat(backup, extension); + + gmx_file_copy(filename, backup, FALSE); + + delete [] backup; + + } + return COLVARS_OK; +} + + +void colvarproxy_gromacs::update_data(const t_commrec *cr, int64_t const step, t_pbc const &pbc, const matrix box, bool bNS) +{ + + if (MASTER(cr)) { + + if(cvm::debug()) { + cvm::log(cvm::line_marker); + cvm::log("colvarproxy_gromacs, step no. "+cvm::to_str(colvars->it)+"\n"+ + "Updating internal data.\n"); + } + + // step update on master only due to the call of colvars pointer. + if (first_timestep) { + first_timestep = false; + } else { + // Use the time step number inherited from GROMACS + if ( step - previous_gmx_step == 1 ) + colvars->it++; + // Other cases? + } + } // end master + + gmx_pbc = pbc; + gmx_box = box; + gmx_bNS = bNS; + + previous_gmx_step = step; + + // Prepare data for MPI communication + if(PAR(cr) && bNS) { + dd_make_local_group_indices(cr->dd->ga2la, n_colvars_atoms, ind, &nat_loc, &ind_loc, &nalloc_loc, xa_ind); + } +} + + +void colvarproxy_gromacs::calculateForces( + const gmx::ForceProviderInput &forceProviderInput, + gmx::ForceProviderOutput *forceProviderOutput) +{ + + const t_commrec *cr = &(forceProviderInput.cr_); + // Local atom coords + const gmx::ArrayRef x = forceProviderInput.x_; + // Local atom coords (coerced into into old gmx type) + const rvec *x_pointer = &(x.data()->as_vec()); + + + // Eventually there needs to be an interface to update local data upon neighbor search + // We could check if by chance all atoms are in one node, and skip communication + communicate_group_positions(cr, x_colvars_unwrapped, xa_shifts, xa_eshifts, + gmx_bNS, x_pointer, n_colvars_atoms, nat_loc, + ind_loc, xa_ind, xa_old_whole, gmx_box); + + // Communicate_group_positions takes care of removing shifts (unwrapping) + // in single node jobs, communicate_group_positions() is efficient and adds no overhead + + if (MASTER(cr)) + { + // On non-master nodes, jump directly to applying the forces + + // Zero the forces on the atoms, so that they can be accumulated by the colvars. + for (size_t i = 0; i < atoms_new_colvar_forces.size(); i++) { + atoms_new_colvar_forces[i].x = atoms_new_colvar_forces[i].y = atoms_new_colvar_forces[i].z = 0.0; + } + + // Get the atom positions from the Gromacs array. + for (size_t i = 0; i < atoms_ids.size(); i++) { + atoms_positions[i] = cvm::rvector(x_colvars_unwrapped[i][0], x_colvars_unwrapped[i][1], x_colvars_unwrapped[i][2]); + } + + bias_energy = 0.0; + // Call the collective variable module to fill atoms_new_colvar_forces + if (colvars->calc() != COLVARS_OK) { + cvm::error("Error calling colvars->calc()\n"); + } + + // Copy the forces to C array for broadcasting + for (int i = 0; i < n_colvars_atoms; i++) + { + f_colvars[i][0] = atoms_new_colvar_forces[i].x; + f_colvars[i][1] = atoms_new_colvar_forces[i].y; + f_colvars[i][2] = atoms_new_colvar_forces[i].z; + } + + forceProviderOutput->enerd_.term[F_COM_PULL] += bias_energy; + } // master node + + //Broadcast the forces to all the nodes + if (PAR(cr)) + { + nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, f_colvars); + } + + const gmx::ArrayRef &f_out = forceProviderOutput->forceWithVirial_.force_; + matrix local_colvars_virial = { { 0 } }; + + // Pass the applied forces back to GROMACS + for (int i = 0; i < n_colvars_atoms; i++) + { + int i_global = ind[i]; + + // check if this is a local atom and find out locndx + if (PAR(cr)) { + const int *locndx = cr->dd->ga2la->findHome(i_global); + if (locndx) { + f_out[*locndx] += f_colvars[i]; + add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); + } + // Do nothing if atom is not local + } else { // Non MPI-parallel + f_out[i_global] += f_colvars[i]; + add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); + } + } + + forceProviderOutput->forceWithVirial_.addVirialContribution(local_colvars_virial); + return; +} + + +void colvarproxy_gromacs::add_virial_term(matrix vir, rvec const f, gmx::RVec const x) +{ + for (int j = 0; j < DIM; j++) { + for (int m = 0; m < DIM; m++) { + vir[j][m] -= 0.5 * f[j] * x[m]; + } + } +} + + +// Pass restraint energy value for current timestep to MD engine +void colvarproxy_gromacs::add_energy (cvm::real energy) +{ + bias_energy += energy; +} + +// **************** ATOMS **************** + +int colvarproxy_gromacs::check_atom_id(int atom_number) +{ + // GROMACS uses zero-based arrays. + int const aid = (atom_number-1); + + if (cvm::debug()) + log("Adding atom "+cvm::to_str(atom_number)+ + " for collective variables calculation.\n"); + + if ( (aid < 0) || (aid >= gmx_atoms.nr) ) { + cvm::error("Error: invalid atom number specified, "+ + cvm::to_str(atom_number)+"\n", COLVARS_INPUT_ERROR); + return COLVARS_INPUT_ERROR; + } + + return aid; +} + + +int colvarproxy_gromacs::init_atom(int atom_number) +{ + // GROMACS uses zero-based arrays. + int aid = atom_number-1; + + for (size_t i = 0; i < atoms_ids.size(); i++) { + if (atoms_ids[i] == aid) { + // this atom id was already recorded + atoms_refcount[i] += 1; + return i; + } + } + + aid = check_atom_id(atom_number); + + if(aid < 0) { + return COLVARS_INPUT_ERROR; + } + + int const index = add_atom_slot(aid); + update_atom_properties(index); + return index; +} + +void colvarproxy_gromacs::update_atom_properties(int index) +{ + + // update mass + double const mass = gmx_atoms.atom[atoms_ids[index]].m; + if (mass <= 0.001) { + this->log("Warning: near-zero mass for atom "+ + cvm::to_str(atoms_ids[index]+1)+ + "; expect unstable dynamics if you apply forces to it.\n"); + } + atoms_masses[index] = mass; + // update charge + atoms_charges[index] = gmx_atoms.atom[atoms_ids[index]].q; +} diff --git a/gromacs/gromacs-2023.x/colvarproxy_gromacs.h b/gromacs/gromacs-2023.x/colvarproxy_gromacs.h new file mode 100644 index 000000000..3f331e82d --- /dev/null +++ b/gromacs/gromacs-2023.x/colvarproxy_gromacs.h @@ -0,0 +1,144 @@ +/// -*- c++ -*- +/* Based on Jeff Comer's old code */ +#ifndef GMX_COLVARS_COLVARPROXY_GROMACS_H +#define GMX_COLVARS_COLVARPROXY_GROMACS_H + +#include "colvarmodule.h" +#include "colvaratoms.h" +#include "colvarproxy.h" +#include "gromacs/random/tabulatednormaldistribution.h" +#include "gromacs/random/threefry.h" +#include "gromacs/mdtypes/forceoutput.h" +#include "gromacs/mdlib/groupcoord.h" +#include "gromacs/mdtypes/iforceprovider.h" +#include "gromacs/mdtypes/observableshistory.h" +#include "gromacs/topology/atoms.h" +#include "colvarproxy_gromacs_version.h" + +/// \brief Communication between colvars and Gromacs (implementation of +/// \link colvarproxy \endlink) +class colvarproxy_gromacs : public colvarproxy, public gmx::IForceProvider { +public: + // GROMACS structures. + //PBC struct + t_pbc gmx_pbc; + //Box + const real (*gmx_box)[3]; + // + t_atoms gmx_atoms; +protected: + bool first_timestep; + + bool colvars_restart; + std::string config_file; + size_t restart_frequency_s; + int previous_gmx_step; + double bias_energy; + + bool gmx_bNS; // Is this a neighbor-search step? Eventually will become unnecessary + + // GROMACS random number generation. + gmx::DefaultRandomEngine rng; // gromacs random number generator + gmx::TabulatedNormalDistribution<> normal_distribution; + + + // Node-local bookkepping data + //! Total number of Colvars atoms + int n_colvars_atoms = 0; + //! Part of the atoms that are local. + int nat_loc = 0; + //! Global indices of the Colvars atoms. + int *ind = nullptr; + //! Local indices of the Colvars atoms. + int *ind_loc = nullptr; + //! Allocation size for ind_loc. + int nalloc_loc = 0; + //! Unwrapped positions for all Colvars atoms, communicated to all nodes. + rvec *x_colvars_unwrapped = nullptr; + //! Shifts for all Colvars atoms, to make molecule(s) whole. + ivec *xa_shifts = nullptr; + //! Extra shifts since last DD step. + ivec *xa_eshifts = nullptr; + //! Old positions for all Colvars atoms on master. + rvec *xa_old_whole = nullptr; + //! Position of each local atom in the collective array. + int *xa_ind = nullptr; + //! Bias forces on all Colvars atoms + rvec *f_colvars = nullptr; +public: + friend class cvm::atom; + colvarproxy_gromacs(); + ~colvarproxy_gromacs(); + + // Initialize colvars. + void init(t_inputrec *gmx_inp, int64_t step, const gmx_mtop_t &mtop, ObservablesHistory* oh, + const std::string &prefix, gmx::ArrayRef filenames_config, + const std::string &filename_restart, const t_commrec *cr, + const rvec x[], rvec **xa_old_whole_colvars_state_p, int *n_colvars_atoms_state_p); + + void dd_make_local_atoms(const t_commrec *cr); + // Called each step before evaluating the force provider + // Should eventually be replaced by the MDmodule interface? + void update_data(const t_commrec *cr, int64_t const step, t_pbc const &pbc, const matrix box, bool bNS); + /*! \brief + * Computes forces. + * + * \param[in] forceProviderInput struct that collects input data for the force providers + * \param[in,out] forceProviderOutput struct that collects output data of the force providers + */ + virtual void calculateForces(const gmx::ForceProviderInput &forceProviderInput, + gmx::ForceProviderOutput *forceProviderOutput); + + // Compute virial tensor for position r and force f, and add to matrix vir + void add_virial_term(matrix vir, rvec const f, gmx::RVec const r); + + void add_energy (cvm::real energy); + void finish(const t_commrec *cr); + + // **************** SYSTEM-WIDE PHYSICAL QUANTITIES **************** + cvm::real rand_gaussian(); + // **************** SIMULATION PARAMETERS **************** + size_t restart_frequency(); + std::string restart_output_prefix(); + std::string output_prefix(); + // **************** PERIODIC BOUNDARY CONDITIONS **************** + cvm::rvector position_distance (cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const; + + // **************** INPUT/OUTPUT **************** + /// Print a message to the main log + void log (std::string const &message); + /// Print a message to the main log and let the rest of the program handle the error + void error (std::string const &message); + /// Print a message to the main log and exit with error code + void fatal_error (std::string const &message); + /// Print a message to the main log and exit normally + void exit (std::string const &message); + /// Request to set the units used internally by Colvars + int set_unit_system(std::string const &units_in, bool colvars_defined); + int backup_file (char const *filename); + /// Read atom identifiers from a file \param filename name of + /// the file (usually a PDB) \param atoms array to which atoms read + /// from "filename" will be appended \param pdb_field (optiona) if + /// "filename" is a PDB file, use this field to determine which are + /// the atoms to be set + int load_atoms (char const *filename, + std::vector &atoms, + std::string const &pdb_field, + double const pdb_field_value = 0.0); + /// Load the coordinates for a group of atoms from a file + /// (usually a PDB); if "pos" is already allocated, the number of its + /// elements must match the number of atoms in "filename" + int load_coords (char const *filename, + std::vector &pos, + const std::vector &indices, + std::string const &pdb_field, + double const pdb_field_value = 0.0); + + int init_atom(int atom_number); + + int check_atom_id(int atom_number); + void update_atom_properties(int index); +}; + +#endif diff --git a/update-colvars-code.sh b/update-colvars-code.sh index 8932391de..affee817e 100755 --- a/update-colvars-code.sh +++ b/update-colvars-code.sh @@ -84,7 +84,7 @@ then elif [ -f "${target}/include/molfile_plugin.h" ] then code="VMD-PLUGINS" -elif [ -f "${target}/src/gromacs/commandline.h" ] +elif [ -f "${target}/src/gromacs/commandline/cmdlineinit.h" ] then code="GROMACS" else @@ -99,7 +99,7 @@ else elif [ -f "${target}/src/VMDApp.h" ] then code="VMD" - elif [ -f "${target}/src/gromacs/commandline.h" ] + elif [ -f "${target}/src/gromacs/commandline/cmdlineinit.h" ] then code="GROMACS" else @@ -197,6 +197,9 @@ then 2022*) GMX_VERSION='2022.x' ;; + 2023*) + GMX_VERSION='2023.x' + ;; *) if [ $force_update = 0 ] ; then echo " ******************************************************************************" @@ -204,6 +207,12 @@ then echo " You may override with -f, but be mindful of compilation or runtime problems." echo " ******************************************************************************" exit 3 + else + echo " ******************************************************************************" + echo " WARNING: Unsupported GROMACS ${GMX_VERSION} version will be patched." + echo " Using GROMACS 2023.x files for patching." + echo " ******************************************************************************" + GMX_VERSION='2023.x' fi ;; esac @@ -342,7 +351,7 @@ then condcopy "${source}/lammps/doc/src/fix_colvars.rst" \ "${target}/doc/src/fix_colvars.rst" rm -f "${target}/doc/src/fix_colvars.txt" - + echo ' done.' if [ ${downloaded_pdf} = 1 ] ; then echo "Note: the PDF manual for the latest Colvars version was downloaded. " @@ -560,7 +569,7 @@ then echo "" # Apply patch for Gromacs files - patch ${patch_opts} -d ${target} < ${source}/gromacs/gromacs-${GMX_VERSION}.patch + # patch ${patch_opts} -d ${target} < ${source}/gromacs/gromacs-${GMX_VERSION}.patch ret_val=$? if [ $ret_val -ne 0 ] then From 510d409d66784490b4be93db4bfbc83a69682473 Mon Sep 17 00:00:00 2001 From: Hubert Date: Fri, 28 Jul 2023 15:48:47 +0200 Subject: [PATCH 2/5] fix typo --- update-colvars-code.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/update-colvars-code.sh b/update-colvars-code.sh index affee817e..d4b2edb96 100755 --- a/update-colvars-code.sh +++ b/update-colvars-code.sh @@ -569,7 +569,7 @@ then echo "" # Apply patch for Gromacs files - # patch ${patch_opts} -d ${target} < ${source}/gromacs/gromacs-${GMX_VERSION}.patch + patch ${patch_opts} -d ${target} < ${source}/gromacs/gromacs-${GMX_VERSION}.patch ret_val=$? if [ $ret_val -ne 0 ] then From 44bf5bc62787c5512e1077c9b5bc0a47d46724ef Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Sun, 20 Aug 2023 11:56:07 +0200 Subject: [PATCH 3/5] More robust way to identify GROMACS source tree --- devel-tools/compile-gromacs.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/devel-tools/compile-gromacs.sh b/devel-tools/compile-gromacs.sh index e20e37cec..788a345ae 100755 --- a/devel-tools/compile-gromacs.sh +++ b/devel-tools/compile-gromacs.sh @@ -21,7 +21,7 @@ compile_gromacs_target() { fi local GMX_SRC_DIR="" - if [ -f "${1}/src/gromacs/commandline.h" ] ; then + if [ -d "${1}/src/gromacs" ] ; then GMX_SRC_DIR=$(realpath "${1}") shift else From be858dc5fb7ca02aa7f465f37c0a0636df6fd84a Mon Sep 17 00:00:00 2001 From: Hubert Date: Mon, 28 Aug 2023 15:55:32 +0200 Subject: [PATCH 4/5] fix Actions for GROMACS 2023 --- .github/workflows/test-backends.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/test-backends.yml b/.github/workflows/test-backends.yml index edaf7caff..325f13a44 100644 --- a/.github/workflows/test-backends.yml +++ b/.github/workflows/test-backends.yml @@ -147,6 +147,7 @@ jobs: backend_name: GROMACS-2023 backend_repo: gromacs/gromacs backend_repo_ref: release-2023 + container_name: CentOS9-devel path_compile_script: devel-tools/compile-gromacs.sh test_lib_directory: gromacs/tests/library - rpath_exe: install/bin/gmx_d \ No newline at end of file + rpath_exe: install/bin/gmx_d From a295d2d3c4012640338dff81071ecbbc7357995c Mon Sep 17 00:00:00 2001 From: Hubert Date: Wed, 30 Aug 2023 09:29:50 +0200 Subject: [PATCH 5/5] fix proxy file --- .../gromacs-2023.x/colvarproxy_gromacs.cpp | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp b/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp index 59b8661f5..6b17cda5a 100644 --- a/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp +++ b/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp @@ -99,12 +99,12 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m // Retrieve the topology of all atoms gmx_atoms = gmx_mtop_global_atoms(mtop); - // Read configuration file and set up the proxy only on the master node. - if (MASTER(cr)) + // Read configuration file and set up the proxy only on the MAIN node. + if (MAIN(cr)) { // initiate module: this object will be the communication proxy - // colvarmodule pointer is only defined on the Master due to the static pointer to colvarproxy. + // colvarmodule pointer is only defined on the MAIN due to the static pointer to colvarproxy. colvars = new colvarmodule(this); version_int = get_version_from_string(COLVARPROXY_VERSION); @@ -140,13 +140,13 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m colvars->it = colvars->it_restart = step; - } // end master + } // end MAIN // MPI initialisation // Initialise attributs for the MPI communication - if(MASTER(cr)) { + if(MAIN(cr)) { // Retrieve the number of colvar atoms n_colvars_atoms = atoms_ids.size(); // Copy their global indices @@ -158,8 +158,8 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m // Let the other nodes know the number of colvar atoms. block_bc(cr->mpi_comm_mygroup, n_colvars_atoms); - // Initialise atoms_new_colvar_forces on non-master nodes - if(!MASTER(cr)) { + // Initialise atoms_new_colvar_forces on non-MAIN nodes + if(!MAIN(cr)) { atoms_new_colvar_forces.reserve(n_colvars_atoms); } } @@ -174,7 +174,7 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m // Prepare data // Manage restart with .cpt - if (MASTER(cr)) + if (MAIN(cr)) { /* colvarsHistory is the struct holding the data saved in the cpt @@ -227,7 +227,7 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m if (PAR(cr)) { nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, xa_old_whole); - snew_bc(MASTER(cr), ind, n_colvars_atoms); + snew_bc(MAIN(cr), ind, n_colvars_atoms); nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, ind); } @@ -245,7 +245,7 @@ void colvarproxy_gromacs::init(t_inputrec *ir, int64_t step, const gmx_mtop_t &m } } - if (MASTER(cr) && cvm::debug()) { + if (MAIN(cr) && cvm::debug()) { cvm::log ("atoms_ids = "+cvm::to_str (atoms_ids)+"\n"); cvm::log ("atoms_refcount = "+cvm::to_str (atoms_refcount)+"\n"); cvm::log ("positions = "+cvm::to_str (atoms_positions)+"\n"); @@ -263,7 +263,7 @@ colvarproxy_gromacs::~colvarproxy_gromacs() void colvarproxy_gromacs::finish(const t_commrec *cr) { - if(MASTER(cr)) { + if(MAIN(cr)) { colvars->write_restart_file(output_prefix_str+".colvars.state"); colvars->write_output_files(); } @@ -381,7 +381,7 @@ int colvarproxy_gromacs::backup_file (char const *filename) void colvarproxy_gromacs::update_data(const t_commrec *cr, int64_t const step, t_pbc const &pbc, const matrix box, bool bNS) { - if (MASTER(cr)) { + if (MAIN(cr)) { if(cvm::debug()) { cvm::log(cvm::line_marker); @@ -389,7 +389,7 @@ void colvarproxy_gromacs::update_data(const t_commrec *cr, int64_t const step, t "Updating internal data.\n"); } - // step update on master only due to the call of colvars pointer. + // step update on MAIN only due to the call of colvars pointer. if (first_timestep) { first_timestep = false; } else { @@ -398,7 +398,7 @@ void colvarproxy_gromacs::update_data(const t_commrec *cr, int64_t const step, t colvars->it++; // Other cases? } - } // end master + } // end MAIN gmx_pbc = pbc; gmx_box = box; @@ -408,7 +408,7 @@ void colvarproxy_gromacs::update_data(const t_commrec *cr, int64_t const step, t // Prepare data for MPI communication if(PAR(cr) && bNS) { - dd_make_local_group_indices(cr->dd->ga2la, n_colvars_atoms, ind, &nat_loc, &ind_loc, &nalloc_loc, xa_ind); + dd_make_local_group_indices(cr->dd->ga2la.get(), n_colvars_atoms, ind, &nat_loc, &ind_loc, &nalloc_loc, xa_ind); } } @@ -434,9 +434,9 @@ void colvarproxy_gromacs::calculateForces( // Communicate_group_positions takes care of removing shifts (unwrapping) // in single node jobs, communicate_group_positions() is efficient and adds no overhead - if (MASTER(cr)) + if (MAIN(cr)) { - // On non-master nodes, jump directly to applying the forces + // On non-MAIN nodes, jump directly to applying the forces // Zero the forces on the atoms, so that they can be accumulated by the colvars. for (size_t i = 0; i < atoms_new_colvar_forces.size(); i++) { @@ -463,7 +463,7 @@ void colvarproxy_gromacs::calculateForces( } forceProviderOutput->enerd_.term[F_COM_PULL] += bias_energy; - } // master node + } // MAIN node //Broadcast the forces to all the nodes if (PAR(cr))