From 1acd9f4b8352f4cef528eddaf0021c4a89aa6eda Mon Sep 17 00:00:00 2001 From: Xu Yang Date: Sat, 26 Sep 2026 19:51:47 +0800 Subject: [PATCH] Feature: PLUMED interface for the MD driver Enable collective variables, biasing potentials and free-energy methods provided by PLUMED (https://www.plumed.org) in ABACUS MD runs: - new INPUT items 'plumed' and 'plumed_file' - the interface is driven from the shared MD loop in run_md.cpp, so all the integrators (nve/nvt/npt/langevin/fire/msst) are supported - ABACUS atomic units are declared to PLUMED through setMD*Units; positions and cell (Bohr), forces (Hartree/Bohr), potential energy (Hartree), masses (amu) are handed over in native units and PLUMED converts them internally - the interface follows the reference call sequence of the PLUMED distribution: values, prepareCalc, setForces/setVirial, performCalc - PLUMED is driven by the master MPI rank only; the biasing forces it adds to the MD force array are broadcast to the other ranks - build with -DPLUMED_ROOT= (./configure --prefix=... && make install) --- CMakeLists.txt | 27 +++ docs/advanced/input_files/input-main.md | 12 ++ source/CMakeLists.txt | 1 + source/module_io/read_input_item_md.cpp | 12 ++ source/module_md/md_base.h | 9 + source/module_md/run_md.cpp | 82 +++++++++ source/module_parameter/md_parameter.h | 3 + source/module_plumed/CMakeLists.txt | 9 + source/module_plumed/plumed_interface.cpp | 206 ++++++++++++++++++++++ source/module_plumed/plumed_interface.h | 118 +++++++++++++ 10 files changed, 479 insertions(+) create mode 100644 source/module_plumed/CMakeLists.txt create mode 100644 source/module_plumed/plumed_interface.cpp create mode 100644 source/module_plumed/plumed_interface.h diff --git a/CMakeLists.txt b/CMakeLists.txt index a2f5dec2bdb..008f5375e92 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -632,6 +632,32 @@ if(DEFINED DeePMD_DIR) endif() endif() +if(DEFINED PLUMED_ROOT) + # enhanced sampling through the PLUMED plugin (https://www.plumed.org); only + # the MD driver makes use of it. PLUMED_ROOT is the prefix of a PLUMED + # installation (created with `./configure --prefix=... && make install`), + # with the library in ${PLUMED_ROOT}/lib and the header in + # ${PLUMED_ROOT}/include/plumed/wrapper/Plumed.h. + message(STATUS "PLUMED support is enabled. PLUMED_ROOT = " ${PLUMED_ROOT}) + add_compile_definitions(USE_PLUMED) + include_directories(${PLUMED_ROOT}/include) + find_library(PLUMED_LIBRARY NAMES plumed PATHS ${PLUMED_ROOT}/lib NO_DEFAULT_PATH) + if(NOT PLUMED_LIBRARY) + message(FATAL_ERROR "PLUMED_ROOT is set but libplumed was not found in ${PLUMED_ROOT}/lib.") + endif() + # The kernel library (a dependency of the wrapper library) is linked + # explicitly so that the executable locates it through its own rpath, + # without requiring LD_LIBRARY_PATH at run time. Its symbols are not + # referenced directly, so it must be protected from the default + # --as-needed of some toolchains. + find_library(PLUMED_KERNEL_LIBRARY NAMES plumedKernel PATHS ${PLUMED_ROOT}/lib NO_DEFAULT_PATH) + target_link_libraries(${ABACUS_BIN_NAME} ${PLUMED_LIBRARY} ${CMAKE_DL_LIBS}) + if(PLUMED_KERNEL_LIBRARY) + target_link_libraries(${ABACUS_BIN_NAME} "-Wl,--no-as-needed" ${PLUMED_KERNEL_LIBRARY} "-Wl,--as-needed") + endif() + set_target_properties(${ABACUS_BIN_NAME} PROPERTIES BUILD_RPATH ${PLUMED_ROOT}/lib) +endif() + if(DEFINED TensorFlow_DIR) find_package(TensorFlow REQUIRED) include_directories(${TensorFlow_DIR}/include) @@ -722,6 +748,7 @@ target_link_libraries( cell symmetry md + plumed_interface planewave surchem neighbor diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 5b4a0232a4c..e8ff0aacfd9 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -2649,6 +2649,18 @@ These variables are used to control molecular dynamics calculations. For more in - **Description**: Whether to output lattice virials into the file `OUT.${suffix}/MD_dump`. - **Default**: True +### plumed + +- **Type**: Boolean +- **Description**: Whether to use PLUMED (https://www.plumed.org) to compute collective variables, biasing potentials and free-energy methods during an MD run. ABACUS must be configured with `-DPLUMED_ROOT=` to enable this feature. +- **Default**: False + +### plumed_file + +- **Type**: String +- **Description**: The input file of PLUMED, read when [`plumed`](#plumed) is enabled. Relative paths are resolved against the working directory of the run. +- **Default**: plumed.dat + ### md_seed - **Type**: Integer diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 1f4d4a83703..a03fd05972f 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -9,6 +9,7 @@ add_subdirectory(module_hsolver) add_subdirectory(module_basis/module_ao) add_subdirectory(module_basis/module_nao) add_subdirectory(module_md) +add_subdirectory(module_plumed) add_subdirectory(module_basis/module_pw) add_subdirectory(module_esolver) add_subdirectory(module_hamilt_lcao/module_gint) diff --git a/source/module_io/read_input_item_md.cpp b/source/module_io/read_input_item_md.cpp index 419a601e282..da9b2a356f0 100644 --- a/source/module_io/read_input_item_md.cpp +++ b/source/module_io/read_input_item_md.cpp @@ -402,5 +402,17 @@ void ReadInput::item_md() read_sync_bool(input.mdp.dump_virial); this->add_item(item); } + { + Input_Item item("plumed"); + item.annotation = "use PLUMED for collective variables and enhanced sampling"; + read_sync_bool(input.mdp.plumed); + this->add_item(item); + } + { + Input_Item item("plumed_file"); + item.annotation = "input file of PLUMED"; + read_sync_string(input.mdp.plumed_file); + this->add_item(item); + } } } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_md/md_base.h b/source/module_md/md_base.h index c48afe10cc4..533a5865c69 100644 --- a/source/module_md/md_base.h +++ b/source/module_md/md_base.h @@ -85,6 +85,15 @@ class MD_base double potential=0.0; ///< potential energy double kinetic; ///< kinetic energy + /** + * @brief get the current time step used by the integrator (in atomic units + * of time, hbar/E_hartree) + */ + double get_md_dt() const + { + return md_dt; + } + protected: const MD_para& mdp; ///< input parameters used in md UnitCell& ucell; ///< unitcell information diff --git a/source/module_md/run_md.cpp b/source/module_md/run_md.cpp index 3b44f627d80..d883d44aaf0 100644 --- a/source/module_md/run_md.cpp +++ b/source/module_md/run_md.cpp @@ -1,5 +1,7 @@ #include "run_md.h" +#include + #include "module_parameter/parameter.h" #include "fire.h" #include "langevin.h" @@ -7,6 +9,7 @@ #include "module_base/global_file.h" #include "module_base/timer.h" #include "module_io/print_info.h" +#include "module_plumed/plumed_interface.h" #include "msst.h" #include "nhchain.h" #include "verlet.h" @@ -14,6 +17,36 @@ namespace Run_MD { +namespace +{ +/// @brief collect the absolute Cartesian positions (in Bohr) of all the atoms +void collect_atomic_positions(const UnitCell& ucell, ModuleBase::Vector3* positions) +{ + int iat = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + /// the UnitCell stores tau in units of the lattice constant lat0 + positions[iat] = ucell.atoms[it].tau[ia] * ucell.lat0; + ++iat; + } + } +} + +/// @brief collect the cell matrix (in Bohr) in the row-major layout expected +/// by PLUMED (row i is the i-th lattice vector) +void collect_cell(const UnitCell& ucell, double* cell) +{ + /// the UnitCell stores the lattice vectors without the lattice constant + const double* latvec = &ucell.latvec.e11; + for (int i = 0; i < 9; ++i) + { + cell[i] = latvec[i] * ucell.lat0; + } +} +} // namespace + void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Parameter& param_in) { ModuleBase::TITLE("Run_MD", "md_line"); @@ -46,12 +79,42 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Paramet ModuleBase::WARNING_QUIT("md_line", "no such md_type!"); } + /// interface to PLUMED: no-op unless `plumed` is set in INPUT and ABACUS + /// was compiled with PLUMED support (-DPLUMED_ROOT=...) + ModulePlumed::PlumedInterface plumed; + + /// absolute Cartesian positions (Bohr) and cell (Bohr) handed over to + /// PLUMED; note that the UnitCell stores tau and latvec in units of the + /// lattice constant lat0 and that mdrun->pos only stores the displacement + /// of the current MD step + std::vector> plumed_pos(unit_in.nat); + double plumed_cell[9]; + /// md cycle while ((mdrun->step_ + mdrun->step_rst_) <= param_in.mdp.md_nstep && !mdrun->stop) { if (mdrun->step_ == 0) { mdrun->setup(p_esolver, PARAM.globalv.global_readin_dir); + + if (param_in.mdp.plumed) + { + plumed.init(param_in.mdp.plumed_file, + unit_in.nat, + mdrun->get_md_dt(), + mdrun->allmass, + param_in.globalv.myrank); + /// the initial configuration already carries positions and + /// forces, so PLUMED can be called right away + collect_atomic_positions(unit_in, plumed_pos.data()); + collect_cell(unit_in, plumed_cell); + plumed.move(mdrun->step_ + mdrun->step_rst_, + plumed_pos.data(), + mdrun->force, + plumed_cell, + mdrun->potential, + param_in.inp.cal_stress ? mdrun->virial.c : nullptr); + } } else { @@ -67,6 +130,21 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Paramet param_in.inp.cal_stress, mdrun->virial); + if (plumed.active()) + { + /// let PLUMED compute the collective variables and add the + /// biasing forces to mdrun->force; they are used by the second + /// half of the current step and by the first half of the next + collect_atomic_positions(unit_in, plumed_pos.data()); + collect_cell(unit_in, plumed_cell); + plumed.move(mdrun->step_ + mdrun->step_rst_, + plumed_pos.data(), + mdrun->force, + plumed_cell, + mdrun->potential, + param_in.inp.cal_stress ? mdrun->virial.c : nullptr); + } + mdrun->second_half(); MD_func::compute_stress(unit_in, @@ -121,7 +199,11 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Paramet mdrun->step_++; } + /// flush the COLVAR and other PLUMED output files + plumed.finalize(); + delete mdrun; + ModuleBase::timer::tick("Run_MD", "md_line"); return; } diff --git a/source/module_parameter/md_parameter.h b/source/module_parameter/md_parameter.h index 80d04b0e94d..9fd1efab1e6 100644 --- a/source/module_parameter/md_parameter.h +++ b/source/module_parameter/md_parameter.h @@ -67,6 +67,9 @@ struct MD_para ///< not. liuyu 2023-03-01 bool dump_virial = true; ///< output lattice virial into the file MD_dump or ///< not. liuyu 2023-03-01 + bool plumed = false; ///< use PLUMED for collective variables and enhanced + ///< sampling (requires ABACUS compiled with PLUMED) + std::string plumed_file = "plumed.dat"; ///< input file of PLUMED }; #endif // MD_PARA_H \ No newline at end of file diff --git a/source/module_plumed/CMakeLists.txt b/source/module_plumed/CMakeLists.txt new file mode 100644 index 00000000000..6b263c85c6d --- /dev/null +++ b/source/module_plumed/CMakeLists.txt @@ -0,0 +1,9 @@ +add_library( + plumed_interface + OBJECT + plumed_interface.cpp +) + +if(ENABLE_COVERAGE) + add_coverage(plumed_interface) +endif() diff --git a/source/module_plumed/plumed_interface.cpp b/source/module_plumed/plumed_interface.cpp new file mode 100644 index 00000000000..5c50c12f674 --- /dev/null +++ b/source/module_plumed/plumed_interface.cpp @@ -0,0 +1,206 @@ +#include "plumed_interface.h" + +#include + +#ifdef __MPI +#include "mpi.h" +#endif + +#include "module_base/tool_quit.h" + +#ifdef USE_PLUMED +// C API of the PLUMED wrapper library, matching in a +// PLUMED installation. The handle is an opaque pointer-sized structure; the +// wrapper library takes care of loading the PLUMED kernel. +extern "C" +{ + struct plumed_handle + { + void* p; + }; + plumed_handle plumed_create(); + void plumed_cmd(plumed_handle plumed, const char* key, const void* val); + void plumed_finalize(plumed_handle plumed); +} +#endif + +namespace ModulePlumed +{ + +namespace +{ +// Conversion factors from ABACUS atomic units to the PLUMED default units +// (length: nm, energy: kJ/mol, time: ps, mass: amu). They are declared to +// PLUMED once in init() with setMD*Units; afterwards all the arrays handed +// over to PLUMED are in ABACUS native units. +constexpr double BOHR_TO_NM = 0.052917721067; // 1 Bohr = 0.052917721067 nm +constexpr double HARTREE_TO_KJMOL = 2625.499638; // 1 Hartree = 2625.499638 kJ/mol +constexpr double AU_TIME_TO_PS = 2.418884326509e-5; // atomic unit of time = 2.418884326509e-5 ps +constexpr double ELECTRON_MASS_TO_AMU = 5.48579909065e-4; // 1 electron mass = 5.48579909065e-4 amu +} // namespace + +PlumedInterface::~PlumedInterface() +{ + this->finalize(); +} + +bool PlumedInterface::init(const std::string& plumed_file, + const int nat, + const double dt_au, + const double* masses_au, + const int my_rank) +{ +#ifdef USE_PLUMED + this->active_run_ = true; + this->nat_ = nat; + this->my_rank_ = my_rank; + + // ABACUS stores atomic masses in atomic units (electron masses), while + // PLUMED expects them in amu + this->masses_amu_.resize(nat); + for (int i = 0; i < nat; ++i) + { + this->masses_amu_[i] = masses_au[i] * ELECTRON_MASS_TO_AMU; + } + + // PLUMED is driven by the master rank only; the biasing forces are + // broadcast to the other ranks in move() + if (my_rank != 0) + { + return true; + } + + this->plumed_ = new plumed_handle(plumed_create()); + if (this->plumed_->p == nullptr) + { + ModuleBase::WARNING_QUIT("ModulePlumed::PlumedInterface::init", + "could not create a PLUMED object, check the PLUMED installation"); + } + + int real_precision = 8; // double precision reals + int natoms = nat; + double length_units = BOHR_TO_NM; + double energy_units = HARTREE_TO_KJMOL; + double time_units = AU_TIME_TO_PS; + double timestep = dt_au; + const char* engine = "abacus"; + + // commands that configure PLUMED must be issued before "init" + plumed_cmd(*this->plumed_, "setRealPrecision", &real_precision); + plumed_cmd(*this->plumed_, "setMDLengthUnits", &length_units); + plumed_cmd(*this->plumed_, "setMDEnergyUnits", &energy_units); + plumed_cmd(*this->plumed_, "setMDTimeUnits", &time_units); + plumed_cmd(*this->plumed_, "setPlumedDat", plumed_file.c_str()); + plumed_cmd(*this->plumed_, "setNatoms", &natoms); + plumed_cmd(*this->plumed_, "setMDEngine", engine); + plumed_cmd(*this->plumed_, "setTimestep", ×tep); + plumed_cmd(*this->plumed_, "init", nullptr); + + return true; +#else + (void)plumed_file; + (void)nat; + (void)dt_au; + (void)masses_au; + (void)my_rank; + + ModuleBase::WARNING_QUIT("ModulePlumed::PlumedInterface::init", + "plumed was requested in INPUT but ABACUS was not compiled with PLUMED support; " + "reconfigure with -DPLUMED_ROOT=/path/to/plumed2"); + return false; +#endif +} + +void PlumedInterface::move(const int istep, + const ModuleBase::Vector3* pos, + ModuleBase::Vector3* force, + const double* cell, + const double potential, + const double* virial) +{ +#ifdef USE_PLUMED + if (!this->active_run_) + { + return; + } + + // Vector3 is a plain {double x, y, z} structure, so the arrays can + // be handed to PLUMED as contiguous [nat][3] buffers + double* force_ptr = reinterpret_cast(force); + + if (this->my_rank_ == 0) + { + double* pos_ptr = reinterpret_cast(const_cast*>(pos)); + int step = istep; + double energy = potential; + double virial_buf[9]; + + plumed_cmd(*this->plumed_, "setStep", &step); + plumed_cmd(*this->plumed_, "setPositions", pos_ptr); + plumed_cmd(*this->plumed_, "setMasses", this->masses_amu_.data()); + plumed_cmd(*this->plumed_, "setBox", cell); + plumed_cmd(*this->plumed_, "setEnergy", &energy); + /// prepare the calculation graph; the arrays that receive the biasing + /// forces are registered afterwards, as done in the reference + /// implementations distributed with PLUMED + plumed_cmd(*this->plumed_, "prepareCalc", nullptr); + plumed_cmd(*this->plumed_, "setForces", force_ptr); + // The Box (cell) data receives forces whenever a periodic collective + // variable is biased, so PLUMED always needs a valid force buffer + // here; when the MD run does not provide a virial, a zero buffer is + // registered instead (as in the reference implementation distributed + // with GROMACS). The virial is handed over in energy units with the + // convention virial_ij = -stress_ij * volume (the same quantity used + // by the official Quantum ESPRESSO interface); ABACUS keeps the + // potential part of the stress (Hartree/Bohr^3) in its MD virial + // matrix, so the volume is multiplied back here. + for (int i = 0; i < 9; ++i) + { + this->virial_buf_[i] = 0.0; + } + if (virial != nullptr) + { + const double volume = std::fabs(cell[0] * (cell[4] * cell[8] - cell[5] * cell[7]) + - cell[1] * (cell[3] * cell[8] - cell[5] * cell[6]) + + cell[2] * (cell[3] * cell[7] - cell[4] * cell[6])); + for (int i = 0; i < 9; ++i) + { + this->virial_buf_[i] = -virial[i] * volume; + } + } + plumed_cmd(*this->plumed_, "setVirial", this->virial_buf_); + plumed_cmd(*this->plumed_, "performCalc", nullptr); + } + +#ifdef __MPI + // PLUMED only ran on the master rank: pass the biasing forces on to all + // the others (the MD force arrays are replicated on every rank) + MPI_Bcast(force_ptr, 3 * this->nat_, MPI_DOUBLE, 0, MPI_COMM_WORLD); +#endif +#else + (void)istep; + (void)pos; + (void)force; + (void)cell; + (void)potential; + (void)virial; +#endif +} + +void PlumedInterface::finalize() +{ +#ifdef USE_PLUMED + if (this->plumed_ != nullptr) + { + if (this->my_rank_ == 0) + { + plumed_finalize(*this->plumed_); + } + delete this->plumed_; + this->plumed_ = nullptr; + } +#endif + this->active_run_ = false; +} + +} // namespace ModulePlumed diff --git a/source/module_plumed/plumed_interface.h b/source/module_plumed/plumed_interface.h new file mode 100644 index 00000000000..bc9214f4002 --- /dev/null +++ b/source/module_plumed/plumed_interface.h @@ -0,0 +1,118 @@ +#ifndef PLUMED_INTERFACE_H +#define PLUMED_INTERFACE_H + +#include +#include + +#include "module_base/vector3.h" + +// opaque handle of the PLUMED wrapper library, see +// in a PLUMED installation +struct plumed_handle; + +namespace ModulePlumed +{ + +/** + * @brief interface between an ABACUS MD run and the PLUMED plugin + * + * This class wraps the C API of the PLUMED wrapper library (``libplumed.so``, + * ``plumed_create`` / ``plumed_cmd`` / ``plumed_finalize``, see + * in a PLUMED installation) so that + * collective variables, biasing potentials and free-energy methods provided + * by PLUMED (https://www.plumed.org) can be used in ABACUS MD runs. + * + * Units: ABACUS atomic units differ from the PLUMED defaults, so the + * conversion factors are declared once in init() through the commands + * ``setMDEnergyUnits`` (Hartree -> kJ/mol), ``setMDLengthUnits`` + * (Bohr -> nm) and ``setMDTimeUnits`` (atomic unit of time -> ps). From then + * on all quantities are passed in ABACUS native units and PLUMED performs the + * conversion internally: + * - positions / cell : Bohr + * - forces : Hartree/Bohr (biasing forces are + * added in place by PLUMED) + * - potential energy : Hartree + * - virial : Hartree (=-stress * volume, + * following the Quantum ESPRESSO + * interface convention) + * - masses : electron masses (converted to amu + * internally) + * + * PLUMED is only driven by the master MPI rank (rank 0); the biasing forces + * are broadcast to all the other ranks at the end of move(). + * + * This interface is a no-op stub unless ABACUS is configured with + * ``-DPLUMED_ROOT=/path/to/plumed2``; requesting PLUMED in INPUT while ABACUS + * was built without it stops the run with an error message. + */ +class PlumedInterface +{ + public: + PlumedInterface() = default; + ~PlumedInterface(); + + PlumedInterface(const PlumedInterface&) = delete; + PlumedInterface& operator=(const PlumedInterface&) = delete; + + /** + * @brief create a PLUMED object and read the PLUMED input file + * + * Must be called once before the MD loop, on all ranks (only rank 0 + * actually creates and initializes PLUMED). + * + * @param plumed_file [in] path of the PLUMED input file + * @param nat [in] number of atoms + * @param dt_au [in] MD time step in atomic units of time + * @param masses_au [in] atomic masses in ABACUS atomic units (electron + * masses), nat entries + * @param my_rank [in] MPI rank of this process + * @return true when the interface is ready to be used in the MD loop + */ + bool init(const std::string& plumed_file, + const int nat, + const double dt_au, + const double* masses_au, + const int my_rank); + + /** + * @brief one MD step: hand over the configuration, let PLUMED compute the + * collective variables and add the biasing forces to ``force`` + * + * @param istep [in] global MD step index + * @param pos [in] atomic positions in Bohr (nat entries) + * @param force [in,out] atomic forces in Hartree/Bohr; PLUMED adds the + * biasing forces in place + * @param cell [in] 3x3 cell matrix in Bohr, row-major, the rows being + * the lattice vectors (or nullptr when not periodic) + * @param potential [in] potential energy in Hartree + * @param virial [in] 3x3 virial in Hartree (=-stress * volume), or + * nullptr when the stress is not available + */ + void move(const int istep, + const ModuleBase::Vector3* pos, + ModuleBase::Vector3* force, + const double* cell, + const double potential, + const double* virial); + + /// @brief finalize PLUMED (flushes the COLVAR and other output files) + void finalize(); + + /// @brief whether PLUMED is active in this run + bool active() const + { + return this->active_run_; + } + + private: + struct plumed_handle* plumed_ = nullptr; ///< handle returned by plumed_create() + bool active_run_ = false; ///< whether init() has been called + int nat_ = 0; ///< number of atoms + int my_rank_ = 0; ///< MPI rank of this process + std::vector masses_amu_; ///< atomic masses in amu + double virial_buf_[9] = {0.0}; ///< virial handed over to PLUMED (energy units) +}; + +} // namespace ModulePlumed + +#endif // PLUMED_INTERFACE_H