From d5041457af0ef865decd7286c083687d24e0f1b9 Mon Sep 17 00:00:00 2001 From: Xu Yang Date: Sun, 27 Sep 2026 04:45:06 +0800 Subject: [PATCH 1/2] Feature: PLUMED interface for the MD driver Port of the PLUMED interface to the MDCell-based MD driver of the develop branch: - new INPUT items 'plumed' and 'plumed_file' - new module source/source_plumed wrapping the PLUMED C API (plumed_create / plumed_cmd / plumed_finalize); the interface is driven from the shared MD loop in source_md/run_md.cpp, so all the integrators are supported - ABACUS atomic units are declared to PLUMED through setMD*Units; the positions/cell (Bohr), forces (Hartree/Bohr), potential energy (Hartree) and masses are handed over in native units, and PLUMED converts them internally - the call sequence follows the reference implementations distributed with PLUMED (values, prepareCalc, setForces/setVirial, performCalc); every step is flagged as a checkpointing step so that PLUMED flushes its output files (COLVAR, ...) while the run is in progress - build with -DPLUMED_ROOT= - current limitation: a single MPI rank is required (the atoms owned by the master rank are handed over), which covers the standard serial/OpenMP MD runs --- CMakeLists.txt | 12 + source/CMakeLists.txt | 20 ++ .../source_io/module_parameter/md_parameter.h | 3 + .../module_parameter/read_input_item_md.cpp | 25 +++ source/source_md/md_base.h | 9 + source/source_md/run_md.cpp | 110 +++++++++ source/source_plumed/CMakeLists.txt | 15 ++ source/source_plumed/plumed_interface.cpp | 212 ++++++++++++++++++ source/source_plumed/plumed_interface.h | 118 ++++++++++ 9 files changed, 524 insertions(+) create mode 100644 source/source_plumed/CMakeLists.txt create mode 100644 source/source_plumed/plumed_interface.cpp create mode 100644 source/source_plumed/plumed_interface.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 33b852f0cd1..0f58c600125 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -617,6 +617,18 @@ if(DEFINED DeePMD_DIR) find_package(DeePMD REQUIRED) endif() +if(DEFINED PLUMED_ROOT) + # PLUMED plugin for collective variables and enhanced sampling + # (https://www.plumed.org). PLUMED_ROOT is the prefix of an installed + # PLUMED (created with `./configure --prefix=... && make install`). + message(STATUS "PLUMED support is enabled. PLUMED_ROOT = " ${PLUMED_ROOT}) + find_library(PLUMED_LIBRARY NAMES plumed PATHS ${PLUMED_ROOT}/lib NO_DEFAULT_PATH) + find_library(PLUMED_KERNEL_LIBRARY NAMES plumedKernel 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() +endif() + if(DEFINED NEP_DIR) find_package(NEP REQUIRED) endif() diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index ffd4272f71f..23c74bff58f 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -292,6 +292,18 @@ if(DEFINED DeePMD_DIR) list(APPEND _abacus_feature_include_dirs ${DeePMD_DIR}/include) endif() +if(DEFINED PLUMED_ROOT) + # PLUMED plugin (https://www.plumed.org). The wrapper library is linked + # together with the kernel library it depends on (protected from the default + # --as-needed of some toolchains so that the kernel also ends up in the + # DT_NEEDED entries of the executable). + list(APPEND _abacus_feature_libs ${PLUMED_LIBRARY} ${CMAKE_DL_LIBS}) + if(PLUMED_KERNEL_LIBRARY) + list(APPEND _abacus_feature_libs "-Wl,--no-as-needed" ${PLUMED_KERNEL_LIBRARY} "-Wl,--as-needed") + endif() + list(APPEND _abacus_feature_include_dirs ${PLUMED_ROOT}/include) +endif() + if(DEFINED NEP_DIR AND NEP_FOUND) list(APPEND _abacus_feature_libs NEP::nep) endif() @@ -538,6 +550,7 @@ add_subdirectory(source_hsolver) add_subdirectory(source_basis/module_ao) add_subdirectory(source_basis/module_nao) add_subdirectory(source_md) +add_subdirectory(source_plumed) add_subdirectory(source_basis/module_pw) add_subdirectory(source_esolver) add_subdirectory(source_io) @@ -688,6 +701,7 @@ target_link_libraries( neighbor_search neighbor md + plumed_interface planewave symmetry cell @@ -739,6 +753,12 @@ if(ENABLE_RAPIDJSON) endif() # External closure: feature consumers precede numerical providers and runtimes. +if(DEFINED PLUMED_ROOT) + # let the executable locate libplumed.so and its kernel library through its + # own rpath (no LD_LIBRARY_PATH needed at run time) + set_target_properties(${ABACUS_BIN_NAME} PROPERTIES BUILD_RPATH ${PLUMED_ROOT}/lib) +endif() + target_link_libraries(${ABACUS_BIN_NAME} PRIVATE abacus::link_libs) # Torch includes some LAPACK routines, but with floating-point exceptions. diff --git a/source/source_io/module_parameter/md_parameter.h b/source/source_io/module_parameter/md_parameter.h index 92a42ce8b65..2bd87bee1f9 100644 --- a/source/source_io/module_parameter/md_parameter.h +++ b/source/source_io/module_parameter/md_parameter.h @@ -71,6 +71,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 diff --git a/source/source_io/module_parameter/read_input_item_md.cpp b/source/source_io/module_parameter/read_input_item_md.cpp index cebe8d7d944..75ab9d0cbe6 100644 --- a/source/source_io/module_parameter/read_input_item_md.cpp +++ b/source/source_io/module_parameter/read_input_item_md.cpp @@ -245,6 +245,31 @@ Note that md_tlast is only used in NVT/NPT simulations. If md_tlast is unset or 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"; + item.category = "Molecular dynamics"; + item.type = "Boolean"; + item.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 " + "compiled with -DPLUMED_ROOT= to enable this feature."; + item.default_value = "False"; + item.unit = ""; + read_sync_bool(input.mdp.plumed); + this->add_item(item); + } + { + Input_Item item("plumed_file"); + item.annotation = "input file of PLUMED"; + item.category = "Molecular dynamics"; + item.type = "String"; + item.description = "The input file of PLUMED, read when `plumed` is enabled. Relative paths are " + "resolved against the working directory of the run."; + item.default_value = "plumed.dat"; + item.unit = ""; + read_sync_string(input.mdp.plumed_file); + this->add_item(item); + } { Input_Item item("md_seed"); item.annotation = "random seed for MD"; diff --git a/source/source_md/md_base.h b/source/source_md/md_base.h index c3f870ee1e8..65a2af84e3a 100644 --- a/source/source_md/md_base.h +++ b/source/source_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 MDCell& mdcell; ///< mdcell information diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index b2eb16420da..71f91f7a8b6 100644 --- a/source/source_md/run_md.cpp +++ b/source/source_md/run_md.cpp @@ -19,8 +19,12 @@ #include "verlet.h" #include "source_cell/update_cell.h" #include "source_cell/print_cell.h" +#include "source_plumed/plumed_interface.h" #include +#ifdef __MPI +#include +#endif namespace Run_MD { @@ -74,6 +78,54 @@ void prepare_mdcell(MDCell& mdcell, UnitCell& ucell, DomainDecomposition& decomp mdcell.mutable_stru_meta() = unitcell::make_stru_meta(ucell); } +namespace +{ +/// @brief collect the absolute Cartesian positions (in Bohr) of the atoms +/// owned by this rank; the cell stores cart in units of the lattice constant +void collect_plumed_positions(const MDCell& mdcell, std::vector>& positions) +{ + const double lat0 = mdcell.lat0(); + const std::vector& atoms = mdcell.owned_atoms(); + positions.resize(atoms.size()); + for (std::size_t i = 0; i < atoms.size(); ++i) + { + positions[i] = atoms[i].cart * lat0; + } +} + +/// @brief gather the atomic forces (in Hartree/Bohr) +void collect_plumed_forces(const MDCell& mdcell, std::vector>& forces) +{ + const std::vector& atoms = mdcell.owned_atoms(); + forces.resize(atoms.size()); + for (std::size_t i = 0; i < atoms.size(); ++i) + { + forces[i] = atoms[i].force; + } +} + +/// @brief write the (PLUMED-biased) forces back into the cell +void apply_plumed_forces(MDCell& mdcell, const std::vector>& forces) +{ + std::vector& atoms = mdcell.owned_atoms(); + for (std::size_t i = 0; i < atoms.size(); ++i) + { + atoms[i].force = forces[i]; + } +} + +/// @brief 3x3 cell matrix (in Bohr, row-major) expected by PLUMED +void collect_plumed_cell(const MDCell& mdcell, double* cell) +{ + const double lat0 = mdcell.lat0(); + const double* latvec = &mdcell.latvec().e11; // unitless lattice vectors + for (int i = 0; i < 9; ++i) + { + cell[i] = latvec[i] * lat0; + } +} +} // namespace + void md_line(MDCell& mdcell, ModuleESolver::ESolver* p_esolver, const Parameter& param_in, @@ -108,12 +160,59 @@ void md_line(MDCell& mdcell, 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; + std::vector> plumed_positions; + std::vector> plumed_forces; + double plumed_cell[9]; + + auto run_plumed_step = [&](const int istep) { + collect_plumed_positions(mdcell, plumed_positions); + collect_plumed_forces(mdcell, plumed_forces); + collect_plumed_cell(mdcell, plumed_cell); + plumed.move(istep, + plumed_positions.data(), + plumed_forces.data(), + plumed_cell, + mdrun->potential, + param_in.inp.cal_stress ? mdrun->virial.c : nullptr); + apply_plumed_forces(mdcell, plumed_forces); + }; + /// md cycle, mohan update 2026-01-04, change '<=' to '<' while ((mdrun->step_ + mdrun->step_rst_) < param_in.mdp.md_nstep && !mdrun->stop) { if (mdrun->step_ == 0) { mdrun->setup(p_esolver, param_in.globalv.global_readin_dir, decomp); + + if (param_in.mdp.plumed) + { + int nproc = 1; +#ifdef __MPI + MPI_Comm_size(MPI_COMM_WORLD, &nproc); +#endif + if (nproc > 1) + { + ModuleBase::WARNING_QUIT("md_line", + "the PLUMED interface currently supports a single MPI rank only; " + "run the MD with one process (OpenMP threads are fine)"); + } + std::vector plumed_masses(mdcell.owned_atoms().size()); + for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) + { + plumed_masses[i] = mdcell.owned_atoms()[i].mass; + } + plumed.init(param_in.mdp.plumed_file, + static_cast(mdcell.owned_atoms().size()), + mdrun->get_md_dt(), + plumed_masses.data(), + param_in.globalv.myrank); + /// the initial configuration already carries positions and + /// forces, so PLUMED can be called right away + run_plumed_step(mdrun->step_ + mdrun->step_rst_); + } } else { @@ -134,6 +233,14 @@ void md_line(MDCell& mdcell, mdrun->virial, param_in.mdp.md_out_force); + if (plumed.active()) + { + /// let PLUMED compute the collective variables and add the + /// biasing forces; they are used by the second half of the + /// current step and by the first half of the next one + run_plumed_step(mdrun->step_ + mdrun->step_rst_); + } + mdrun->second_half(); MD_func::compute_stress(mdcell, @@ -172,6 +279,9 @@ void md_line(MDCell& mdcell, mdrun->step_++; } + /// flush the COLVAR and other PLUMED output files + plumed.finalize(); + delete mdrun; ModuleBase::timer::end("Run_MD", "md_line"); return; diff --git a/source/source_plumed/CMakeLists.txt b/source/source_plumed/CMakeLists.txt new file mode 100644 index 00000000000..7def1dfbb93 --- /dev/null +++ b/source/source_plumed/CMakeLists.txt @@ -0,0 +1,15 @@ +add_library( + plumed_interface + OBJECT + plumed_interface.cpp +) + +# The PLUMED interface is only compiled in when the build is configured with +# -DPLUMED_ROOT= +if(DEFINED PLUMED_ROOT) + target_compile_definitions(plumed_interface PRIVATE USE_PLUMED) +endif() + +if(ENABLE_COVERAGE) + add_coverage(plumed_interface) +endif() diff --git a/source/source_plumed/plumed_interface.cpp b/source/source_plumed/plumed_interface.cpp new file mode 100644 index 00000000000..c2ae9e04f1e --- /dev/null +++ b/source/source_plumed/plumed_interface.cpp @@ -0,0 +1,212 @@ +#include "plumed_interface.h" + +#include + +#ifdef __MPI +#include "mpi.h" +#endif + +#include "source_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_); + // flag this step as a checkpointing step: PLUMED then flushes all its + // output files (COLVAR, ...) in update(), so that the collective + // variables can be followed while the MD is running. Without it, + // PLUMED only flushes every 10000 steps. + int check_point = 1; + plumed_cmd(*this->plumed_, "doCheckPoint", &check_point); + 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/source_plumed/plumed_interface.h b/source/source_plumed/plumed_interface.h new file mode 100644 index 00000000000..9f8d8165a45 --- /dev/null +++ b/source/source_plumed/plumed_interface.h @@ -0,0 +1,118 @@ +#ifndef PLUMED_INTERFACE_H +#define PLUMED_INTERFACE_H + +#include +#include + +#include "source_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 From b4fd38916be6011dfdb28d09ba26a85dc596a8e8 Mon Sep 17 00:00:00 2001 From: Xu Yang Date: Wed, 30 Sep 2026 12:01:17 +0800 Subject: [PATCH 2/2] docs: document the PLUMED interface --- docs/advanced/input_files/input-main.md | 14 +++++ docs/advanced/interface/index.rst | 1 + docs/advanced/interface/plumed.md | 69 +++++++++++++++++++++++++ docs/parameters.yaml | 16 ++++++ 4 files changed, 100 insertions(+) create mode 100644 docs/advanced/interface/plumed.md diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 428f38fd01d..c87b57a7852 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -363,6 +363,8 @@ - [dump\_force](#dump_force) - [dump\_vel](#dump_vel) - [dump\_virial](#dump_virial) + - [plumed](#plumed) + - [plumed\_file](#plumed_file) - [md\_seed](#md_seed) - [md\_tfreq](#md_tfreq) - [md\_tchain](#md_tchain) @@ -3598,6 +3600,18 @@ - **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 compiled with -DPLUMED_ROOT=<prefix of a PLUMED installation> to enable this feature. +- **Default**: False + +### plumed_file + +- **Type**: String +- **Description**: The input file of PLUMED, read when `plumed` is enabled. Relative paths are resolved against the working directory of the run. +- **Default**: plumed.dat + ### md_seed - **Type**: Integer diff --git a/docs/advanced/interface/index.rst b/docs/advanced/interface/index.rst index 6622e821b7e..1671fcd671e 100644 --- a/docs/advanced/interface/index.rst +++ b/docs/advanced/interface/index.rst @@ -19,3 +19,4 @@ Interfaces to Other Softwares ShengBTE candela TB2J + plumed diff --git a/docs/advanced/interface/plumed.md b/docs/advanced/interface/plumed.md new file mode 100644 index 00000000000..5b5d3eacea9 --- /dev/null +++ b/docs/advanced/interface/plumed.md @@ -0,0 +1,69 @@ +# PLUMED + +## Introduction + +[PLUMED](https://www.plumed.org) is an open-source library for molecular +dynamics that provides collective variables (CVs), biasing potentials +(metadynamics, umbrella sampling, walls, ...) and free-energy analysis tools. +It is driven through a small C API, so the same `plumed.dat` input file can be +reused across different MD engines. + +ABACUS supports PLUMED during *ab initio* molecular dynamics (`calculation md`): +on every MD step the atomic positions, cell, masses, potential energy and (when +available) the virial are handed over to PLUMED, and the bias forces returned +by PLUMED are added to the nuclear forces. The interface is driven from the +shared MD loop, so all the algorithms available through `md_type` (`nve`, +`nvt`, ...) are supported. + +> **Note:** the current implementation requires a **single MPI rank** (OpenMP +> parallelization within that rank is fully supported). Running PLUMED together +> with multi-rank MPI is not supported yet. + +## Installation + +Install PLUMED 2.x following the instructions on https://www.plumed.org/doc +(`./configure --prefix= && make install`), then configure +ABACUS with the prefix of that installation: + +```bash +cmake -B build -DPLUMED_ROOT= +cmake --build build -j +``` + +When `PLUMED_ROOT` is not given, the feature is disabled and nothing else +changes, so PLUMED is a completely optional dependency. + +## How to use + +Enable PLUMED in the `INPUT` file and point to the PLUMED input file: + +``` +plumed 1 +plumed_file plumed.dat +``` + +A minimal `plumed.dat` for well-tempered metadynamics along a distance: + +``` +UNITS LENGTH=A TIME=fs ENERGY=eV +d: DISTANCE ATOMS=1,2 +metad: METAD ARG=d PACE=100 HEIGHT=0.05 SIGMA=0.1 BIASFACTOR=10 TEMP=300 FILE=HILLS +PRINT ARG=d,metad.bias FILE=COLVAR STRIDE=10 +``` + +The output files requested by PLUMED (`COLVAR`, `HILLS`, ...) are written into +the working directory of the run and are flushed at every MD step, so they can +be followed while the simulation is running. The accumulated bias can be +converted into a free-energy surface with the PLUMED command line tool, e.g. + +```bash +plumed sum_hills --hills HILLS --mintozero --outfile fes.dat +``` + +### Units + +The interface declares the ABACUS atomic units (Bohr, Hartree and the atomic +unit of time) to PLUMED, and all the arrays are exchanged in these units; +PLUMED converts them internally, so no manual conversion is needed. Use the +`UNITS` directive in `plumed.dat` to select the units of the values written in +the file (for example `UNITS LENGTH=A TIME=fs ENERGY=eV`). diff --git a/docs/parameters.yaml b/docs/parameters.yaml index a0f76236e86..cd0d6bff18f 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -1536,6 +1536,22 @@ parameters: default_value: "True" unit: "" availability: "" + - name: plumed + category: Molecular dynamics + 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 compiled with -DPLUMED_ROOT= to enable this feature. + default_value: "False" + unit: "" + availability: "" + - name: plumed_file + category: Molecular dynamics + type: String + description: | + The input file of PLUMED, read when `plumed` is enabled. Relative paths are resolved against the working directory of the run. + default_value: plumed.dat + unit: "" + availability: "" - name: md_seed category: Molecular dynamics type: Integer