Feature: PLUMED interface for the MD driver

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()

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

docs/advanced/interface/index.rst

@@ -19,3 +19,4 @@ Interfaces to Other Softwares
    ShengBTE
    candela
    TB2J
+   plumed

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=<plumed-install> && make install`), then configure
+ABACUS with the prefix of that installation:
+
+```bash
+cmake -B build -DPLUMED_ROOT=<plumed-install>
+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`).

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=<prefix of a PLUMED installation> 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

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.

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

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=<prefix of a PLUMED installation> 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";

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