From c2ca37f46ef9d1e71ad5c9280b4f4ed50541ba66 Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Fri, 4 Oct 2024 09:01:38 -0400 Subject: [PATCH] Remove patch files for old and unsupported GROMACS versions --- .github/workflows/test-backends.yml | 12 - gromacs/gromacs-2018.x.patch | 453 -------------- gromacs/gromacs-2020.x.patch | 501 --------------- .../gromacs-2020.x/colvarproxy_gromacs.cpp | 578 ----------------- gromacs/gromacs-2020.x/colvarproxy_gromacs.h | 144 ----- gromacs/gromacs-2021.x.patch | 509 --------------- .../gromacs-2021.x/colvarproxy_gromacs.cpp | 580 ------------------ gromacs/gromacs-2021.x/colvarproxy_gromacs.h | 144 ----- gromacs/gromacs-2022.x.patch | 521 ---------------- gromacs/gromacs-2022.x/CMakeLists.txt | 35 -- .../gromacs-2022.x/colvarproxy_gromacs.cpp | 580 ------------------ gromacs/gromacs-2022.x/colvarproxy_gromacs.h | 144 ----- gromacs/gromacs-2023.x.patch | 513 ---------------- gromacs/gromacs-2023.x/CMakeLists.txt | 35 -- .../gromacs-2023.x/colvarproxy_gromacs.cpp | 580 ------------------ gromacs/gromacs-2023.x/colvarproxy_gromacs.h | 144 ----- 16 files changed, 5473 deletions(-) delete mode 100644 gromacs/gromacs-2018.x.patch delete mode 100644 gromacs/gromacs-2020.x.patch delete mode 100644 gromacs/gromacs-2020.x/colvarproxy_gromacs.cpp delete mode 100644 gromacs/gromacs-2020.x/colvarproxy_gromacs.h delete mode 100644 gromacs/gromacs-2021.x.patch delete mode 100644 gromacs/gromacs-2021.x/colvarproxy_gromacs.cpp delete mode 100644 gromacs/gromacs-2021.x/colvarproxy_gromacs.h delete mode 100644 gromacs/gromacs-2022.x.patch delete mode 100644 gromacs/gromacs-2022.x/CMakeLists.txt delete mode 100644 gromacs/gromacs-2022.x/colvarproxy_gromacs.cpp delete mode 100644 gromacs/gromacs-2022.x/colvarproxy_gromacs.h delete mode 100644 gromacs/gromacs-2023.x.patch delete mode 100644 gromacs/gromacs-2023.x/CMakeLists.txt delete mode 100644 gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp delete mode 100644 gromacs/gromacs-2023.x/colvarproxy_gromacs.h diff --git a/.github/workflows/test-backends.yml b/.github/workflows/test-backends.yml index fd311ffd1..c946855e1 100644 --- a/.github/workflows/test-backends.yml +++ b/.github/workflows/test-backends.yml @@ -68,18 +68,6 @@ jobs: private_key: ${{ secrets.PULL_VMD_KEY }} private_key_vmd_plugins: ${{ secrets.PULL_VMD_PLUGINS_KEY }} - gromacs-2023: - name: GROMACS (release-2023) - uses: ./.github/workflows/backend-template.yml - with: - 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_mpi_d - gromacs-main: name: GROMACS (main) uses: ./.github/workflows/backend-template.yml diff --git a/gromacs/gromacs-2018.x.patch b/gromacs/gromacs-2018.x.patch deleted file mode 100644 index 068425a5f..000000000 --- a/gromacs/gromacs-2018.x.patch +++ /dev/null @@ -1,453 +0,0 @@ -diff -aur src/gromacs/CMakeLists.txt src_colvars/gromacs/CMakeLists.txt ---- src/gromacs/CMakeLists.txt 2020-10-08 11:47:51.703335032 +0200 -+++ src_colvars/gromacs/CMakeLists.txt 2020-10-08 11:51:21.327529496 +0200 -@@ -105,6 +105,7 @@ - add_subdirectory(awh) - add_subdirectory(simd) - add_subdirectory(imd) -+add_subdirectory(colvars) - if (NOT GMX_BUILD_MDRUN_ONLY) - add_subdirectory(gmxana) - add_subdirectory(gmxpreprocess) -Seulement dans src_colvars/gromacs: colvars -diff -aur src/gromacs/fileio/checkpoint.cpp src_colvars/gromacs/fileio/checkpoint.cpp ---- src/gromacs/fileio/checkpoint.cpp 2020-10-08 11:47:51.751335992 +0200 -+++ src_colvars/gromacs/fileio/checkpoint.cpp 2020-10-08 11:51:21.295528856 +0200 -@@ -72,6 +72,7 @@ - #include "gromacs/mdtypes/observableshistory.h" - #include "gromacs/mdtypes/state.h" - #include "gromacs/mdtypes/swaphistory.h" -+#include "gromacs/mdtypes/colvarshistory.h" - #include "gromacs/trajectory/trajectoryframe.h" - #include "gromacs/utility/arrayref.h" - #include "gromacs/utility/baseversion.h" -@@ -99,7 +100,10 @@ - * But old code can not read a new entry that is present in the file - * (but can read a new format when new entries are not present). - */ --static const int cpt_version = 17; -+/* -+ * Add a value of 1000 for Colvars support -+ */ -+static const int cpt_version = 17 + 1000; - - - const char *est_names[estNR] = -@@ -886,7 +887,7 @@ - int *natoms, int *ngtc, int *nnhpres, int *nhchainlength, - int *nlambda, int *flags_state, - int *flags_eks, int *flags_enh, int *flags_dfh, int *flags_awhh, -- int *nED, int *eSwapCoords, -+ int *nED, int *eSwapCoords, int *ecolvars, - FILE *list) - { - bool_t res = 0; -@@ -1046,6 +1047,14 @@ - { - *flags_awhh = 0; - } -+ if (*file_version >= 17 + 1000) -+ { -+ do_cpt_int_err(xd, "colvars atoms", ecolvars, list); -+ } -+ else -+ { -+ *ecolvars = 0; -+ } - } - - static int do_cpt_footer(XDR *xd, int file_version) -@@ -1512,6 +1521,37 @@ - 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, - gmx::CorrelationGridHistory *corrGrid, - FILE *list, int eawhh) -@@ -1918,13 +1958,17 @@ - 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); -+ - do_cpt_header(gmx_fio_getxdr(fp), FALSE, &file_version, - &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime, - &eIntegrator, &simulation_part, &step, &t, &nppnodes, - DOMAINDECOMP(cr) ? domdecCells : nullptr, &npmenodes, - &state->natoms, &state->ngtc, &state->nnhpres, - &state->nhchainlength, &nlambda, &state->flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, -- &nED, &eSwapCoords, -+ &nED, &eSwapCoords, &ecolvars, - nullptr); - - sfree(version); -@@ -1940,6 +1984,7 @@ - (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, nED, edsamhist, nullptr) < 0) || - (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), NULL) < 0) || - (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, eSwapCoords, swaphist, nullptr) < 0) || -+ (do_cpt_colvars(gmx_fio_getxdr(fp), FALSE, ecolvars, colvarshist, nullptr) < 0) || - (do_cpt_files(gmx_fio_getxdr(fp), FALSE, &outputfiles, &noutputfiles, nullptr, - file_version) < 0)) - { -@@ -2197,7 +2242,7 @@ - int eIntegrator_f, nppnodes_f, npmenodes_f; - ivec dd_nc_f; - int natoms, ngtc, nnhpres, nhchainlength, nlambda, fflags, flags_eks, flags_enh, flags_dfh, flags_awhh; -- int nED, eSwapCoords; -+ int nED, eSwapCoords, ecolvars; - int ret; - gmx_file_position_t *outputfiles; - int nfiles; -@@ -2228,7 +2273,7 @@ - &nppnodes_f, dd_nc_f, &npmenodes_f, - &natoms, &ngtc, &nnhpres, &nhchainlength, &nlambda, - &fflags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, -- &nED, &eSwapCoords, nullptr); -+ &nED, &eSwapCoords, &ecolvars, nullptr); - - if (bAppendOutputFiles && - file_version >= 13 && double_prec != GMX_DOUBLE) -@@ -2415,6 +2460,16 @@ - cp_error(); - } - -+ if (ecolvars != 0 && observablesHistory->colvarsHistory == nullptr) -+ { -+ observablesHistory->colvarsHistory = std::unique_ptr(new colvarshistory_t{}); -+ } -+ ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, ecolvars, observablesHistory->colvarsHistory.get(), nullptr); -+ if (ret) -+ { -+ cp_error(); -+ } -+ - ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, &nfiles, nullptr, file_version); - if (ret) - { -@@ -2630,7 +2685,7 @@ - int flags_eks, flags_enh, flags_dfh, flags_awhh; - double t; - t_state state; -- int nED, eSwapCoords; -+ int nED, eSwapCoords, ecolvars; - t_fileio *fp; - - if (filename == nullptr || -@@ -2651,7 +2706,7 @@ - &eIntegrator, simulation_part, step, &t, &nppnodes, dd_nc, &npme, - &state.natoms, &state.ngtc, &state.nnhpres, &state.nhchainlength, - &nlambda, &state.flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, -- &nED, &eSwapCoords, nullptr); -+ &nED, &eSwapCoords, &ecolvars, nullptr); - - gmx_fio_close(fp); - } -@@ -2668,7 +2723,7 @@ - ivec dd_nc; - int nlambda; - int flags_eks, flags_enh, flags_dfh, flags_awhh; -- int nED, eSwapCoords; -+ int nED, eSwapCoords, ecolvars; - int nfiles_loc; - gmx_file_position_t *files_loc = nullptr; - int ret; -@@ -2678,7 +2733,7 @@ - &eIntegrator, simulation_part, step, t, &nppnodes, dd_nc, &npme, - &state->natoms, &state->ngtc, &state->nnhpres, &state->nhchainlength, - &nlambda, &state->flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, -- &nED, &eSwapCoords, nullptr); -+ &nED, &eSwapCoords, &ecolvars, nullptr); - ret = - do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr); - if (ret) -@@ -2725,6 +2780,13 @@ - cp_error(); - } - -+ colvarshistory_t colvarshist = {}; -+ ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, ecolvars, &colvarshist, nullptr); -+ if (ret) -+ { -+ cp_error(); -+ } -+ - ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, - &files_loc, - &nfiles_loc, -@@ -2823,7 +2885,7 @@ - ivec dd_nc; - int nlambda; - int flags_eks, flags_enh, flags_dfh, flags_awhh;; -- int nED, eSwapCoords; -+ int nED, eSwapCoords, ecolvars; - int ret; - gmx_file_position_t *outputfiles; - int nfiles; -@@ -2836,7 +2898,7 @@ - &eIntegrator, &simulation_part, &step, &t, &nppnodes, dd_nc, &npme, - &state.natoms, &state.ngtc, &state.nnhpres, &state.nhchainlength, - &nlambda, &state.flags, -- &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, &nED, &eSwapCoords, -+ &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, &nED, &eSwapCoords, &ecolvars, - out); - ret = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out); - if (ret) -@@ -2878,6 +2940,12 @@ - } - - if (ret == 0) -+ { -+ colvarshistory_t colvarshist = {}; -+ ret = do_cpt_colvars(gmx_fio_getxdr(fp), TRUE, ecolvars, &colvarshist, out); -+ } -+ -+ if (ret == 0) - { - ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, &nfiles, out, file_version); - } -diff -aur src/gromacs/mdlib/sim_util.cpp src_colvars/gromacs/mdlib/sim_util.cpp ---- src/gromacs/mdlib/sim_util.cpp 2020-10-08 11:47:51.775336472 +0200 -+++ src_colvars/gromacs/mdlib/sim_util.cpp 2020-10-08 11:51:21.307529096 +0200 -@@ -111,6 +111,7 @@ - #include "gromacs/utility/pleasecite.h" - #include "gromacs/utility/smalloc.h" - #include "gromacs/utility/sysinfo.h" -+#include "gromacs/colvars/colvarproxy_gromacs.h" - - #include "nbnxn_gpu.h" - #include "nbnxn_kernels/nbnxn_kernel_cpu.h" -@@ -839,10 +840,22 @@ - */ - if (computeForces) - { -+ - /* Collect forces from modules */ - forceProviders->calculateForces(cr, mdatoms, box, t, x, forceWithVirial); - } - -+ /* COLVARS */ -+ /* Colvars Module needs some updated data then compute the forces */ -+ if (inputrec->bColvars) -+ { -+ t_pbc pbc; -+ set_pbc(&pbc, inputrec->ePBC, box); -+ inputrec->colvars_proxy->update_data(cr, step, pbc, box, bNS); -+ -+ inputrec->colvars_proxy->calculateForces(cr, box, x, forceWithVirial); -+ } -+ - if (inputrec->bPull && pull_have_potential(inputrec->pull_work)) - { - pull_potential_wrapper(cr, inputrec, box, x, -diff -aur src/gromacs/mdtypes/colvarshistory.h src_colvars/gromacs/mdtypes/colvarshistory.h ---- src/gromacs/mdtypes/colvarshistory.h 1970-01-01 01:00:00.000000000 +0100 -+++ src_colvars/gromacs/mdtypes/colvarshistory.h 2020-10-08 11:51:21.131525574 +0200 -@@ -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 -aur src/gromacs/mdtypes/inputrec.h src_colvars/gromacs/mdtypes/inputrec.h ---- src/gromacs/mdtypes/inputrec.h 2020-10-08 11:47:51.571332390 +0200 -+++ src_colvars/gromacs/mdtypes/inputrec.h 2020-10-08 11:51:21.131525574 +0200 -@@ -50,6 +50,7 @@ - #define EGP_TABLE (1<<1) - - struct pull_params_t; -+class colvarproxy_gromacs; - - namespace gmx - { -@@ -398,6 +399,10 @@ - gmx_bool useTwinRange; // Whether twin-range scheme is active - always false if a valid .tpr was read - - gmx::KeyValueTreeObject *params; -+ -+ /* COLVARS */ -+ gmx_bool bColvars; /* Do we do colvars calculations ? */ -+ colvarproxy_gromacs *colvars_proxy; /* The object for the colvars calculations */ - }; - - int ir_optimal_nstcalcenergy(const t_inputrec *ir); -diff -aur src/gromacs/mdtypes/observableshistory.cpp src_colvars/gromacs/mdtypes/observableshistory.cpp ---- src/gromacs/mdtypes/observableshistory.cpp 2020-10-08 11:47:51.571332390 +0200 -+++ src_colvars/gromacs/mdtypes/observableshistory.cpp 2020-10-08 11:51:21.131525574 +0200 -@@ -40,6 +40,7 @@ - #include "gromacs/mdtypes/edsamhistory.h" - #include "gromacs/mdtypes/energyhistory.h" - #include "gromacs/mdtypes/swaphistory.h" -+#include "gromacs/mdtypes/colvarshistory.h" - - ObservablesHistory::ObservablesHistory() = default; - ObservablesHistory::~ObservablesHistory() = default; -diff -aur src/gromacs/mdtypes/observableshistory.h src_colvars/gromacs/mdtypes/observableshistory.h ---- src/gromacs/mdtypes/observableshistory.h 2020-10-08 11:47:51.571332390 +0200 -+++ src_colvars/gromacs/mdtypes/observableshistory.h 2020-10-08 11:51:21.131525574 +0200 -@@ -58,6 +58,7 @@ - class energyhistory_t; - struct edsamhistory_t; - struct swaphistory_t; -+struct colvarshistory_t; - - /*! \libinternal \brief Observables history, for writing/reading to/from checkpoint file - */ -@@ -72,6 +73,9 @@ - //! Ion/water position swapping history - std::unique_ptr swapHistory; - -+ //! Colvars -+ std::unique_ptr colvarsHistory; -+ - //! Default constructor - ObservablesHistory(); - -diff -aur src/programs/mdrun/md.cpp src_colvars/programs/mdrun/md.cpp ---- src/programs/mdrun/md.cpp 2020-10-08 11:47:51.507331110 +0200 -+++ src_colvars/programs/mdrun/md.cpp 2020-10-08 11:51:21.071524374 +0200 -@@ -126,6 +126,9 @@ - #include "gromacs/utility/real.h" - #include "gromacs/utility/smalloc.h" - -+#include "gromacs/colvars/colvarproxy_gromacs.h" -+#include "gromacs/mdtypes/colvarshistory.h" -+ - #include "deform.h" - #include "membed.h" - #include "repl_ex.h" -@@ -837,6 +840,52 @@ - fprintf(fplog, "\n"); - } - -+ /* COLVARS */ -+ if (opt2bSet("-colvars",nfile,fnm)) -+ { -+ -+ char **fnms; -+ int nfile_in; -+ std::string filename_restart; -+ std::string prefix; -+ -+ ir->bColvars = TRUE; -+ -+ /* Retrieve filenames */ -+ nfile_in = opt2fns(&fnms, "-colvars", nfile, fnm); -+ std::vector< std::string> filenames(fnms, fnms+nfile_in); -+ -+ if (opt2bSet("-colvars_restart",nfile,fnm)) -+ { -+ filename_restart = opt2fn("-colvars_restart",nfile,fnm); -+ } -+ -+ /* Determine the prefix for the colvars output files, based on the logfile name. */ -+ std::string logfile = ftp2fn(efLOG, nfile, fnm); -+ /* 4 = ".log".length() */ -+ if(logfile.length() > 4) -+ { -+ prefix = logfile.substr(0,logfile.length()-4); -+ } -+ -+ ir->colvars_proxy = new colvarproxy_gromacs(); -+ ir->colvars_proxy->init(ir,ir->init_step,top_global, observablesHistory, prefix, filenames,filename_restart, cr, MASTER(cr) ? as_rvec_array(state->x.data()) : nullptr); -+ } -+ else -+ { -+ ir->bColvars = FALSE; -+ if (opt2bSet("-colvars_restart",nfile,fnm)) -+ { -+ 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."); -+ } -+ } -+ - walltime_accounting_start(walltime_accounting); - wallcycle_start(wcycle, ewcRUN); - print_start(fplog, cr, walltime_accounting, "mdrun"); -@@ -2025,6 +2074,13 @@ - delete ir->awh; - } - -+ /* COLVARS */ -+ if (ir->bColvars) -+ { -+ ir->colvars_proxy->finish(cr); -+ delete ir->colvars_proxy; -+ } -+ - // Clean up swapcoords - if (ir->eSwapCoords != eswapNO) - { -diff -aur src/programs/mdrun/runner.h src_colvars/programs/mdrun/runner.h ---- src/programs/mdrun/runner.h 2020-10-08 11:47:51.499330950 +0200 -+++ src_colvars/programs/mdrun/runner.h 2020-10-08 11:51:21.071524374 +0200 -@@ -87,7 +87,7 @@ - //! Parallelism-related user options. - gmx_hw_opt_t hw_opt; - //! Filenames and properties from command-line argument values. -- std::array filenames = -+ std::array filenames = - {{{ efTPR, nullptr, nullptr, ffREAD }, - { efTRN, "-o", nullptr, ffWRITE }, - { efCOMPRESSED, "-x", nullptr, ffOPTWR }, -@@ -121,6 +121,8 @@ - { efTOP, "-mp", "membed", ffOPTRD }, - { efNDX, "-mn", "membed", ffOPTRD }, - { efXVG, "-if", "imdforces", ffOPTWR }, -+ { efDAT, "-colvars", "colvars", ffOPTRDMULT }, /* COLVARS */ -+ { efDAT, "-colvars_restart", "colvars", ffOPTRD }, /* COLVARS */ - { efXVG, "-swap", "swapions", ffOPTWR }}}; - /*! \brief Filename arguments. - * diff --git a/gromacs/gromacs-2020.x.patch b/gromacs/gromacs-2020.x.patch deleted file mode 100644 index baef003be..000000000 --- a/gromacs/gromacs-2020.x.patch +++ /dev/null @@ -1,501 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -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) -+ - # The earliest version of the CUDA toolkit that supports c++14 is 9.0 - 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 ---- 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() -+list(APPEND libgromacs_object_library_dependencies colvars) -+list(APPEND libgromacs_object_library_dependencies lepton) -+ - configure_file(version.h.cmakein version.h) - if(GMX_INSTALL_LEGACY_API) - install(FILES -@@ -195,6 +201,8 @@ else() - add_library(libgromacs ${LIBGROMACS_SOURCES}) - endif() - -+gmx_include_colvars_headers() -+ - # Add these contents first because linking their tests can take a lot - # of time, so we want lots of parallel work still available after - # linking starts. -diff --git a/src/gromacs/fileio/checkpoint.cpp b/src/gromacs/fileio/checkpoint.cpp -index 9c6cfe4213..ec2e64f61f 100644 ---- a/src/gromacs/fileio/checkpoint.cpp -+++ b/src/gromacs/fileio/checkpoint.cpp -@@ -73,6 +73,7 @@ - #include "gromacs/mdtypes/pullhistory.h" - #include "gromacs/mdtypes/state.h" - #include "gromacs/mdtypes/swaphistory.h" -+#include "gromacs/mdtypes/colvarshistory.h" - #include "gromacs/trajectory/trajectoryframe.h" - #include "gromacs/utility/arrayref.h" - #include "gromacs/utility/baseversion.h" -@@ -127,7 +128,9 @@ enum cptv - * Backward compatibility for reading old run input files is maintained - * by checking this version number against that of the file and then using - * the correct code path. */ --static const int cpt_version = cptv_Count - 1; -+/* 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 - { - contents->flagsPullHistory = 0; - } -+ if (contents->file_version >= cptv_Count - 1 + 1000) -+ { -+ do_cpt_int_err(xd, "colvars atoms", &contents->ecolvars, list); -+ } -+ else -+ { -+ contents->ecolvars = 0; -+ } -+ - } - - 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) -+{ -+ -+ 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_whole"); -+ 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, -@@ -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); -+ - CheckpointHeaderContents headerContents = { 0, - { 0 }, - { 0 }, -@@ -2357,7 +2402,8 @@ void write_checkpoint(const char* fn, - flags_dfh, - flags_awhh, - nED, -- eSwapCoords }; -+ eSwapCoords, -+ ecolvars }; - std::strcpy(headerContents.version, gmx_version()); - std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath()); - std::strcpy(headerContents.ftime, timebuf.c_str()); -@@ -2377,6 +2423,7 @@ void write_checkpoint(const char* fn, - || (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, nED, edsamhist, nullptr) < 0) - || (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), nullptr) < 0) - || (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, eSwapCoords, swaphist, nullptr) < 0) -+ || (do_cpt_colvars(gmx_fio_getxdr(fp), FALSE, ecolvars, colvarshist, 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?"); -@@ -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{}); -+ } -+ 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) -@@ -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) -+ { -+ cp_error(); -+ } -+ - 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 = {}; -+ 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 fb8f7268be..6feb181b30 100644 ---- a/src/gromacs/fileio/checkpoint.h -+++ b/src/gromacs/fileio/checkpoint.h -@@ -175,6 +175,8 @@ struct CheckpointHeaderContents - int nED; - //! Enum for coordinate swapping. - int eSwapCoords; -+ //! 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 ---- a/src/gromacs/mdlib/energyoutput.cpp -+++ b/src/gromacs/mdlib/energyoutput.cpp -@@ -238,7 +238,7 @@ EnergyOutput::EnergyOutput(ener_file* fp_ene, - bEner_[F_DISPCORR] = (ir->eDispCorr != edispcNO); - 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] = ((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 -index f2528d78b4..7f91cac83d 100644 ---- a/src/gromacs/mdlib/sim_util.cpp -+++ b/src/gromacs/mdlib/sim_util.cpp -@@ -114,6 +114,8 @@ - #include "gromacs/utility/strconvert.h" - #include "gromacs/utility/sysinfo.h" - -+#include "colvarproxy_gromacs.h" -+ - using gmx::AtomLocality; - using gmx::DomainLifetimeWorkload; - using gmx::ForceOutputs; -@@ -553,6 +555,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->ePBC, box); -+ inputrec->colvars_proxy->update_data(cr, step, pbc, box, didNeighborSearch); -+ } -+ - 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 -+++ b/src/gromacs/mdrun/legacymdrunoptions.h -@@ -121,7 +121,9 @@ 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 9ff4b3817d..eb31f1fa89 100644 ---- a/src/gromacs/mdrun/replicaexchange.cpp -+++ b/src/gromacs/mdrun/replicaexchange.cpp -@@ -611,6 +611,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 c2b3c088d7..bb38d44ab2 100644 ---- a/src/gromacs/mdrun/runner.cpp -+++ b/src/gromacs/mdrun/runner.cpp -@@ -115,6 +115,7 @@ - #include "gromacs/mdtypes/md_enums.h" - #include "gromacs/mdtypes/mdrunoptions.h" - #include "gromacs/mdtypes/observableshistory.h" -+#include "gromacs/mdtypes/colvarshistory.h" - #include "gromacs/mdtypes/simulation_workload.h" - #include "gromacs/mdtypes/state.h" - #include "gromacs/mdtypes/state_propagator_data_gpu.h" -@@ -156,6 +157,8 @@ - #include "gromacs/utility/smalloc.h" - #include "gromacs/utility/stringutil.h" - -+#include "colvarproxy_gromacs.h" -+ - #include "isimulator.h" - #include "replicaexchange.h" - #include "simulatorbuilder.h" -@@ -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())) -+ { -+ -+ 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,inputrec->init_step,&mtop, &observablesHistory, prefix, filenames_colvars,filename_restart, cr, MASTER(cr) ? globalState->x.rvec_array() : nullptr, -+ MASTER(cr) ? &globalState->xa_old_whole_colvars : nullptr, MASTER(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 (DOMAINDECOMP(cr)) - { - GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP"); -@@ -1654,6 +1702,16 @@ int Mdrunner::mdrunner() - free_membed(membed); - } - -+ /* 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 a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h -new file mode 100644 -index 0000000000..ea69e03419 ---- /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 copy of xa_old_whole in colvarproxy_gromacs -+} colvarshistory_t; -+ -+#endif -diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 266670f3ef..f1d287d35d 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; -+ -+ /* COLVARS */ -+ 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 ---- a/src/gromacs/mdtypes/observableshistory.cpp -+++ b/src/gromacs/mdtypes/observableshistory.cpp -@@ -41,6 +41,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 d2ba1d820f..a5747139d7 100644 ---- a/src/gromacs/mdtypes/observableshistory.h -+++ b/src/gromacs/mdtypes/observableshistory.h -@@ -59,6 +59,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 - */ -@@ -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 ---- a/src/gromacs/mdtypes/state.cpp -+++ b/src/gromacs/mdtypes/state.cpp -@@ -240,7 +240,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) - - { - // 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 -index a54bff29bb..8619c7935a 100644 ---- a/src/gromacs/mdtypes/state.h -+++ 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 ---- 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-2020.x/colvarproxy_gromacs.cpp b/gromacs/gromacs-2020.x/colvarproxy_gromacs.cpp deleted file mode 100644 index a2203f562..000000000 --- a/gromacs/gromacs-2020.x/colvarproxy_gromacs.cpp +++ /dev/null @@ -1,578 +0,0 @@ -/// -*- c++ -*- - -#include -#include -#include -#include - - -#include "gromacs/math/units.h" -#include "gromacs/mdlib/mdatoms.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,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_ = 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; - // Don't strip the input_prefix_str because colvarmodule.cpp doesn't know that restart file from GROMACS needs the .dat extension. - } - - // 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, 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(); - - /* 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, n_colvars_atoms, xa_old_whole); - snew_bc(cr, ind, n_colvars_atoms); - nblock_bc(cr, 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 { - 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, n_colvars_atoms, f_colvars); - } - - const gmx::ArrayRef &f_out = forceProviderOutput->forceWithVirial_.force_; - matrix local_colvars_virial = { { 0 } }; - const bool computeVirial = forceProviderOutput->forceWithVirial_.computeVirial_; - - // 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]; - if (computeVirial) { - 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]; - if (computeVirial) { - add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); - } - } - } - - if (computeVirial) { - 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-2020.x/colvarproxy_gromacs.h b/gromacs/gromacs-2020.x/colvarproxy_gromacs.h deleted file mode 100644 index 186dd0fa3..000000000 --- a/gromacs/gromacs-2020.x/colvarproxy_gromacs.h +++ /dev/null @@ -1,144 +0,0 @@ -/// -*- 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; - int64_t 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, 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/gromacs/gromacs-2021.x.patch b/gromacs/gromacs-2021.x.patch deleted file mode 100644 index 3f8bdf9f9..000000000 --- a/gromacs/gromacs-2021.x.patch +++ /dev/null @@ -1,509 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index 3715ce1656..65134d0568 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -572,6 +572,9 @@ include(gmxManageTNG) - - include(gmxManageLmfit) - -+include(gmxManageColvars) -+include(gmxManageLepton) -+ - # The earliest version of the CUDA toolkit that supports c++14 is 9.0 - 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 ---- 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() -+list(APPEND libgromacs_object_library_dependencies colvars) -+list(APPEND libgromacs_object_library_dependencies lepton) -+ - configure_file(version.h.cmakein version.h) - if(GMX_INSTALL_LEGACY_API) - install(FILES -@@ -189,6 +195,8 @@ else() - add_library(libgromacs ${LIBGROMACS_SOURCES}) - endif() - -+gmx_include_colvars_headers() -+ - # Add these contents first because linking their tests can take a lot - # of time, so we want lots of parallel work still available after - # linking starts. -diff --git a/src/gromacs/fileio/checkpoint.cpp b/src/gromacs/fileio/checkpoint.cpp -index fadc3d283e..8a6ce56a6f 100644 ---- a/src/gromacs/fileio/checkpoint.cpp -+++ b/src/gromacs/fileio/checkpoint.cpp -@@ -72,6 +72,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" -@@ -177,7 +178,9 @@ enum cptv - * Backward compatibility for reading old run input files is maintained - * by checking this version number against that of the file and then using - * the correct code path. */ --static const int cpt_version = cptv_Count - 1; -+/* 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 - { - contents->isModularSimulatorCheckpoint = false; - } -+ if (contents->file_version >= cptv_Count - 1 + 1000) -+ { -+ do_cpt_int_err(xd, "colvars atoms", &contents->ecolvars, list); -+ } -+ else -+ { -+ contents->ecolvars = 0; -+ } -+ - } - - 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) -+{ -+ -+ 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, -@@ -2330,6 +2371,7 @@ void write_checkpoint_data(t_fileio* fp, - || (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, headerContents.eSwapCoords, - 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?"); -@@ -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{}); -+ } -+ 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) -@@ -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) -+ { -+ cp_error(); -+ } -+ - 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 = {}; -+ 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 c3dbc3c107..c32dfa4f44 100644 ---- a/src/gromacs/fileio/checkpoint.h -+++ b/src/gromacs/fileio/checkpoint.h -@@ -238,6 +238,8 @@ struct CheckpointHeaderContents - int nED; - //! Enum for coordinate swapping. - int 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 0b72fc4e0c..60666542b4 100644 ---- a/src/gromacs/mdlib/energyoutput.cpp -+++ b/src/gromacs/mdlib/energyoutput.cpp -@@ -242,7 +242,7 @@ EnergyOutput::EnergyOutput(ener_file* fp_ene, - bEner_[F_DISPCORR] = (ir->eDispCorr != edispcNO); - 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] = ((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 -index 7e06e4fc9e..38f5604ebd 100644 ---- a/src/gromacs/mdlib/mdoutf.cpp -+++ b/src/gromacs/mdlib/mdoutf.cpp -@@ -64,6 +64,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" -@@ -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); -+ - CheckpointHeaderContents headerContents = { 0, - { 0 }, - { 0 }, -@@ -385,6 +390,7 @@ static void write_checkpoint(const char* fn, - 0, - nED, - eSwapCoords, -+ ecolvars, - false }; - 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 2571b0d216..4572e2baea 100644 ---- a/src/gromacs/mdlib/sim_util.cpp -+++ b/src/gromacs/mdlib/sim_util.cpp -@@ -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, - */ - 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, 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 -+++ b/src/gromacs/mdrun/legacymdrunoptions.h -@@ -127,7 +127,9 @@ 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 c40161d9ef..490db3f10f 100644 ---- a/src/gromacs/mdrun/replicaexchange.cpp -+++ b/src/gromacs/mdrun/replicaexchange.cpp -@@ -610,6 +610,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 232d994e1a..8937b83296 100644 ---- a/src/gromacs/mdrun/runner.cpp -+++ b/src/gromacs/mdrun/runner.cpp -@@ -125,6 +125,7 @@ - #include "gromacs/mdtypes/mdatom.h" - #include "gromacs/mdtypes/mdrunoptions.h" - #include "gromacs/mdtypes/observableshistory.h" -+#include "gromacs/mdtypes/colvarshistory.h" - #include "gromacs/mdtypes/simulation_workload.h" - #include "gromacs/mdtypes/state.h" - #include "gromacs/mdtypes/state_propagator_data_gpu.h" -@@ -167,6 +168,8 @@ - #include "gromacs/utility/smalloc.h" - #include "gromacs/utility/stringutil.h" - -+#include "colvarproxy_gromacs.h" -+ - #include "isimulator.h" - #include "membedholder.h" - #include "replicaexchange.h" -@@ -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())) -+ { -+ -+ 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, MASTER(cr) ? globalState->x.rvec_array() : nullptr, -+ MASTER(cr) ? &globalState->xa_old_whole_colvars : nullptr, MASTER(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 (DOMAINDECOMP(cr)) - { - GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP"); -@@ -1839,6 +1887,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 a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h -new file mode 100644 -index 0000000000..6605e6fce2 ---- /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/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 6e4ee727ab..042acdc01b 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; -+ -+ /* COLVARS */ -+ 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 ---- a/src/gromacs/mdtypes/observableshistory.cpp -+++ b/src/gromacs/mdtypes/observableshistory.cpp -@@ -41,6 +41,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 d2ba1d820f..a5747139d7 100644 ---- a/src/gromacs/mdtypes/observableshistory.h -+++ b/src/gromacs/mdtypes/observableshistory.h -@@ -59,6 +59,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 - */ -@@ -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 ---- a/src/gromacs/mdtypes/state.cpp -+++ b/src/gromacs/mdtypes/state.cpp -@@ -302,7 +302,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) - - { - // 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 -index e38f3f7dbc..06a1cd8484 100644 ---- a/src/gromacs/mdtypes/state.h -+++ 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 ---- 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-2021.x/colvarproxy_gromacs.cpp b/gromacs/gromacs-2021.x/colvarproxy_gromacs.cpp deleted file mode 100644 index 0dfaac12a..000000000 --- a/gromacs/gromacs-2021.x/colvarproxy_gromacs.cpp +++ /dev/null @@ -1,580 +0,0 @@ -/// -*- 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,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_ = 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; - // Don't strip the input_prefix_str because colvarmodule.cpp doesn't know that restart file from GROMACS needs the .dat extension. - } - - // 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 { - 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 } }; - const bool computeVirial = forceProviderOutput->forceWithVirial_.computeVirial_; - - // 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]; - if (computeVirial) { - 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]; - if (computeVirial) { - add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); - } - } - } - - if (computeVirial) { - 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-2021.x/colvarproxy_gromacs.h b/gromacs/gromacs-2021.x/colvarproxy_gromacs.h deleted file mode 100644 index 186dd0fa3..000000000 --- a/gromacs/gromacs-2021.x/colvarproxy_gromacs.h +++ /dev/null @@ -1,144 +0,0 @@ -/// -*- 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; - int64_t 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, 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/gromacs/gromacs-2022.x.patch b/gromacs/gromacs-2022.x.patch deleted file mode 100644 index c4cf77ab2..000000000 --- a/gromacs/gromacs-2022.x.patch +++ /dev/null @@ -1,521 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index cd748c9..31e3b62 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -588,6 +588,10 @@ include(gmxManageLmfit) - - include(gmxManageMuparser) - -+include(gmxManageColvars) -+ -+include(gmxManageLepton) -+ - ################################################## - # Process SIMD instruction settings - ################################################## -diff --git a/src/gromacs/CMakeLists.txt b/src/gromacs/CMakeLists.txt -index 99cc804..40fbe29 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) -+list(APPEND libgromacs_object_library_dependencies lepton) -+ - if(GMX_INSTALL_LEGACY_API) - install(FILES - analysisdata.h -@@ -184,6 +189,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 34114cd..b469541 100644 ---- a/src/gromacs/applied_forces/CMakeLists.txt -+++ b/src/gromacs/applied_forces/CMakeLists.txt -@@ -73,6 +73,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 f6764a1..279953d 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" -@@ -1287,6 +1293,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) -@@ -2035,6 +2049,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, -@@ -2432,6 +2476,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?"); -@@ -2825,6 +2870,17 @@ static void read_checkpoint(const char* 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) -@@ -3000,6 +3056,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) -@@ -3120,6 +3183,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 = {}; -+ 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 5145e1d..3054462 100644 ---- a/src/gromacs/fileio/checkpoint.h -+++ b/src/gromacs/fileio/checkpoint.h -@@ -220,6 +220,8 @@ enum class CheckPointVersion : int - MDModules, - //! Added checkpointing for modular simulator. - ModularSimulator, -+ //! 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 6c0fca5..e06e511 100644 ---- a/src/gromacs/mdlib/energyoutput.cpp -+++ b/src/gromacs/mdlib/energyoutput.cpp -@@ -254,7 +254,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 0a8653a..b75c8b1 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 }, -@@ -381,6 +386,7 @@ static void write_checkpoint(const char* fn, - 0, - nED, - eSwapCoords, -+ ecolvars, - false }; - 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 ---- 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, - */ - 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..b37d822 100644 ---- a/src/gromacs/mdrun/legacymdrunoptions.h -+++ b/src/gromacs/mdrun/legacymdrunoptions.h -@@ -124,7 +124,9 @@ 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 9fa553c..9c5be82 100644 ---- a/src/gromacs/mdrun/replicaexchange.cpp -+++ b/src/gromacs/mdrun/replicaexchange.cpp -@@ -617,6 +617,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 ea7f402..d3ffd59 100644 ---- a/src/gromacs/mdrun/runner.cpp -+++ b/src/gromacs/mdrun/runner.cpp -@@ -126,6 +126,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" -@@ -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" - #include "membedholder.h" - #include "replicaexchange.h" -@@ -2055,6 +2058,52 @@ 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, MASTER(cr) ? globalState->x.rvec_array() : nullptr, -+ MASTER(cr) ? &globalState->xa_old_whole_colvars : nullptr, MASTER(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"); -@@ -2225,6 +2274,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 a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h -new file mode 100644 -index 0000000000..6605e6fce2 ---- /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/inputrec.h b/src/gromacs/mdtypes/inputrec.h -index 922a22e..9d981ff 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; -+ -+ /* COLVARS */ -+ 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 ---- 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 99b315e..d49a386 100644 ---- a/src/gromacs/mdtypes/state.cpp -+++ b/src/gromacs/mdtypes/state.cpp -@@ -369,7 +369,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) - - { - // 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 -index 579e343..725e841 100644 ---- a/src/gromacs/mdtypes/state.h -+++ 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 ---- 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-2022.x/CMakeLists.txt b/gromacs/gromacs-2022.x/CMakeLists.txt deleted file mode 100644 index 551c252d3..000000000 --- a/gromacs/gromacs-2022.x/CMakeLists.txt +++ /dev/null @@ -1,35 +0,0 @@ -# -# 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-2022.x/colvarproxy_gromacs.cpp b/gromacs/gromacs-2022.x/colvarproxy_gromacs.cpp deleted file mode 100644 index 07cc387a9..000000000 --- a/gromacs/gromacs-2022.x/colvarproxy_gromacs.cpp +++ /dev/null @@ -1,580 +0,0 @@ -/// -*- 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; - // Don't strip the input_prefix_str because colvarmodule.cpp doesn't know that restart file from GROMACS needs the .dat extension. - } - - // 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 { - 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 } }; - const bool computeVirial = forceProviderOutput->forceWithVirial_.computeVirial_; - - // 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]; - if (computeVirial) { - 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]; - if (computeVirial) { - add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); - } - } - } - - if (computeVirial) { - 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-2022.x/colvarproxy_gromacs.h b/gromacs/gromacs-2022.x/colvarproxy_gromacs.h deleted file mode 100644 index 4aeb7bc5e..000000000 --- a/gromacs/gromacs-2022.x/colvarproxy_gromacs.h +++ /dev/null @@ -1,144 +0,0 @@ -/// -*- 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; - int64_t 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/gromacs/gromacs-2023.x.patch b/gromacs/gromacs-2023.x.patch deleted file mode 100644 index 2b1cfecb2..000000000 --- a/gromacs/gromacs-2023.x.patch +++ /dev/null @@ -1,513 +0,0 @@ -diff --git a/CMakeLists.txt b/CMakeLists.txt -index 294258d59d..dac9691eff 100644 ---- a/CMakeLists.txt -+++ b/CMakeLists.txt -@@ -638,6 +638,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 8ddc9c19ad..21b32212c9 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 114bfd9834..ba8a3e8dbd 100644 ---- a/src/gromacs/CMakeLists.txt -+++ b/src/gromacs/CMakeLists.txt -@@ -411,6 +411,11 @@ gmx_manage_lmfit() - target_link_libraries(libgromacs PRIVATE lmfit) - target_link_libraries(libgromacs PRIVATE muparser::muparser) - -+gmx_manage_colvars() -+target_link_libraries(libgromacs PRIVATE colvars) -+gmx_manage_lepton() -+target_link_libraries(libgromacs PRIVATE lepton) -+ - # Make sure we fix "everything" found by compilers that support that - gmx_warn_on_everything(libgromacs) - if (CMAKE_CXX_COMPILER_ID STREQUAL "MSVC") -diff --git a/src/gromacs/applied_forces/CMakeLists.txt b/src/gromacs/applied_forces/CMakeLists.txt -index 3c4987f892..53e6b91d70 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 29bb86bf07..fb279dd0ea 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 4811cb883a..d2a32d01d4 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 6fbd4238b0..c57b4dbefe 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 3ea5b7322c..1a2a51be9f 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 }, -@@ -381,6 +386,7 @@ static void write_checkpoint(const char* fn, - 0, - nED, - eSwapCoords, -+ ecolvars, - false }; - std::strcpy(headerContents.version, gmx_version()); - std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath().u8string().c_str()); -diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp -index 5ef0f97f96..112db78878 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 2b6079950e..05606a2cdf 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 6da922e0ef..89088f3c1d 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 387b5fff47..e18423ef8a 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 a/src/gromacs/mdtypes/colvarshistory.h b/src/gromacs/mdtypes/colvarshistory.h -new file mode 100644 -index 0000000000..6605e6fce2 ---- /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 12230eed9c..13f5df1030 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 6ed0e3fc00..172ead081f 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 6a38311e1d..872673f624 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 c8c05e83f0..2026f26186 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 d86e4871c9..9285e7d6be 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 deleted file mode 100644 index 551c252d3..000000000 --- a/gromacs/gromacs-2023.x/CMakeLists.txt +++ /dev/null @@ -1,35 +0,0 @@ -# -# 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 deleted file mode 100644 index 3a3cd1d08..000000000 --- a/gromacs/gromacs-2023.x/colvarproxy_gromacs.cpp +++ /dev/null @@ -1,580 +0,0 @@ -/// -*- 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; - // Don't strip the input_prefix_str because colvarmodule.cpp doesn't know that restart file from GROMACS needs the .dat extension. - } - - // 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 MAIN node. - if (MAIN(cr)) - { - - // initiate module: this object will be the communication proxy - // 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); - - 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 MAIN - - - // MPI initialisation - - // Initialise attributs for the MPI communication - if(MAIN(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-MAIN nodes - if(!MAIN(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 (MAIN(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(MAIN(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 (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"); - 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(MAIN(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 (MAIN(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 MAIN only due to the call of colvars pointer. - if (first_timestep) { - first_timestep = false; - } else { - if ( step - previous_gmx_step == 1 ) - colvars->it++; - // Other cases? - } - } // end MAIN - - 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.get(), 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 (MAIN(cr)) - { - // 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++) { - 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; - } // MAIN 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 } }; - const bool computeVirial = forceProviderOutput->forceWithVirial_.computeVirial_; - - // 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]; - if (computeVirial) { - 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]; - if (computeVirial) { - add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]); - } - } - } - - if (computeVirial) { - 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 deleted file mode 100644 index 4aeb7bc5e..000000000 --- a/gromacs/gromacs-2023.x/colvarproxy_gromacs.h +++ /dev/null @@ -1,144 +0,0 @@ -/// -*- 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; - int64_t 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