[Feature] Add NEP as esolver (Useful information to know that NEP has been interfaced with ABACUS)

CMakeLists.txt

@@ -654,6 +654,15 @@ if(DEFINED DeePMD_DIR)
   endif()
 endif()
 
+if(DEFINED NEP_DIR)
+  find_package(NEP REQUIRED)
+
+  if(NEP_FOUND)
+    add_compile_definitions(__NEP)
+    target_link_libraries(${ABACUS_BIN_NAME} NEP::nep_cpu)
+  endif()
+endif()
+
 if(DEFINED TensorFlow_DIR)
   find_package(TensorFlow REQUIRED)
   include_directories(${TensorFlow_DIR}/include)

cmake/FindNEP.cmake

@@ -0,0 +1,47 @@
+###############################################################################
+# - Find NEP_CPU
+# Finds the NEP_CPU header and library.
+#
+# This module will search for the NEP_CPU library, looking for a hint
+# from the NEP_DIR environment variable or CMake variable.
+#
+# This module defines the following variables:
+#
+#  NEP_FOUND        - True if the NEP_CPU library and headers were found.
+#  NEP_INCLUDE_DIR  - The directory where nep.h is located.
+#  NEP_LIBRARY      - The full path to the NEP_CPU library.
+#
+# It also defines the following imported target:
+#
+#  NEP::nep_cpu     - The NEP_CPU library target.
+#
+###############################################################################
+
+find_path(NEP_INCLUDE_DIR nep.h
+    HINTS ${NEP_DIR}
+    PATH_SUFFIXES "include"
+    )
+
+find_library(NEP_LIBRARY
+    NAMES nepcpu
+    HINTS ${NEP_DIR}
+    PATH_SUFFIXES "lib"
+    )
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(NEP
+    DEFAULT_MSG
+    NEP_LIBRARY NEP_INCLUDE_DIR)
+
+if(NEP_FOUND)
+    if(NOT TARGET NEP::nep_cpu)
+        add_library(NEP::nep_cpu UNKNOWN IMPORTED)
+        set_target_properties(NEP::nep_cpu PROPERTIES
+            IMPORTED_LINK_INTERFACE_LANGUAGES "C"
+            IMPORTED_LOCATION "${NEP_LIBRARY}"
+            INTERFACE_INCLUDE_DIRECTORIES "${NEP_INCLUDE_DIR}"
+        )
+    endif()
+endif()
+
+mark_as_advanced(NEP_INCLUDE_DIR NEP_LIBRARY)

docs/advanced/input_files/input-main.md

@@ -535,6 +535,7 @@ These variables are used to control general system parameters.
   - tddft: real-time time-dependent density functional theory (TDDFT)
   - lj: Leonard Jones potential
   - dp: DeeP potential, see details in [md.md](../md.md#dpmd)
+  - nep: Neuroevolution Potential, see details in [md.md](../md.md#nep)
   - ks-lr: Kohn-Sham density functional theory + LR-TDDFT (Under Development Feature)
   - lr: LR-TDDFT with given KS orbitals (Under Development Feature)
 - **Default**: ksdft

docs/advanced/install.md

@@ -55,6 +55,17 @@ Similarly, DeePMD-kit supports PyTorch backend but its libraries are placed at a
 cmake -B build -DDeePMD_DIR=/dir_to_deepmd-kit -DTorch_DIR=/dir_to_pytorch
 ```
 
+## Build with NEP
+This interface enables running MD simulations with the NEP model. It requires the [NEP_CPU](https://github.com/brucefan1983/NEP_CPU) library, which can be easily installed using toolchain as shown below:
+```bash
+./install_abacus_toolchain.sh --with-nepcpu=install
+```
+
+To build ABACUS:
+```bash
+cmake -B build -DNEP_DIR=/path/to/nep_cpu
+```
+
 ## Build with LibRI and LibComm
 
 The new EXX implementation depends on two external libraries:

docs/advanced/md.md

@@ -87,4 +87,8 @@ ABACUS performs the [Multi-Scale Shock Technique (MSST) integration](https://jou
 Compiling ABACUS with [DeePMD-kit](https://github.com/deepmodeling/deepmd-kit), MD calculations based on machine learning DP model is enabled.
 
 To employ DPMD calculations, [esolver_type](./input_files/input-main.md#esolver_type) should be set to `dp`.
-And the filename of DP model is specified by keyword [pot_file](./input_files/input-main.md#pot_file).
+And the filename of DP model is specified by keyword [pot_file](./input_files/input-main.md#pot_file).
+
+## NEP
+
+If ABACUS is compiled with the Neuroevolution Potential ([NEP](https://gpumd.org/potentials/nep.html)), MD simulations using  NEP models are enabled. To use this feature, set [esolver_type](./input_files/input-main.md#esolver_type) to `nep` and specify the potential file path with the [pot_file](./input_files/input-main.md#pot_file) keyword in your INPUT file.

source/Makefile.Objects

@@ -267,6 +267,7 @@ OBJS_ESOLVER=esolver.o\
     esolver_sdft_pw.o\
     esolver_lj.o\
     esolver_dp.o\
+    esolver_nep.o\
     esolver_of.o\
     esolver_of_tddft.o\
     esolver_of_tool.o\

source/source_esolver/CMakeLists.txt

@@ -7,6 +7,7 @@ list(APPEND objects
     esolver_sdft_pw.cpp
     esolver_lj.cpp
     esolver_dp.cpp
+    esolver_nep.cpp
     esolver_of.cpp
     esolver_of_tddft.cpp
     esolver_of_interface.cpp

source/source_esolver/esolver.cpp

@@ -18,6 +18,7 @@ extern "C"
 }
 #endif
 #include "esolver_dp.h"
+#include "esolver_nep.h"
 #include "esolver_lj.h"
 #include "esolver_of.h"
 #include "esolver_of_tddft.h"
@@ -97,6 +98,10 @@ std::string determine_type()
     {
         esolver_type = "dp_pot";
     }
+    else if (PARAM.inp.esolver_type == "nep")
+    {
+        esolver_type = "nep_pot";
+    }
     else if (esolver_type == "none")
     {
         ModuleBase::WARNING_QUIT("ESolver", "No such esolver_type combined with basis_type");
@@ -338,6 +343,10 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
     {
         return new ESolver_DP(PARAM.mdp.pot_file);
     }
+    else if (esolver_type == "nep_pot")
+    {
+        return new ESolver_NEP(PARAM.mdp.pot_file);
+    }
     throw std::invalid_argument("esolver_type = " + std::string(esolver_type) + ". Wrong in " + std::string(__FILE__)
                                 + " line " + std::to_string(__LINE__));
 }

source/source_esolver/esolver_nep.cpp

@@ -0,0 +1,213 @@
+/**
+ * @file esolver_nep.cpp
+#include "source_io/module_parameter/parameter.h"
+ * @brief Implementation of ESolver_NEP class for neuroevolution potential (NEP).
+ *
+ * This file contains the implementation of the ESolver_NEP class, which is used for solving the energy and forces in a
+ * NEP simulation.
+ * NEP is a method for training deep neural networks to accurately predict the potential energy surface of a
+ * molecular system.
+ *
+ * For more information about NEP, see the following reference:
+ * 1. https://gpumd.org/potentials/nep.html
+ * 2. https://doi.org/10.1002/mgea.70028
+ *
+ * @author MoseyQAQ
+ * @date 2025-10-10
+ */
+#include "esolver_nep.h"
+
+#include "source_base/parallel_common.h"
+#include "source_base/timer.h"
+#include "source_io/output_log.h"
+#include "source_io/cif_io.h"
+
+#include <numeric>
+#include <unordered_map>
+
+using namespace ModuleESolver;
+
+void ESolver_NEP::before_all_runners(UnitCell& ucell, const Input_para& inp)
+{   
+    nep_potential = 0.0;
+    nep_force.create(ucell.nat, 3);
+    nep_virial.create(3, 3);
+    atype.resize(ucell.nat);
+    _e.resize(ucell.nat);
+    _f.resize(3 * ucell.nat);
+    _v.resize(9 * ucell.nat);
+
+    ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", 
+                               ucell, 
+                               "# Generated by ABACUS ModuleIO::CifParser",
+                               "data_?");
+
+#ifdef __NEP
+    /// determine the type map from STRU to NEP model
+    type_map(ucell);
+#endif
+}
+
+void ESolver_NEP::runner(UnitCell& ucell, const int istep)
+{
+    ModuleBase::TITLE("ESolver_NEP", "runner");
+    ModuleBase::timer::tick("ESolver_NEP", "runner");
+
+    // note that NEP are column major, thus a transpose is needed
+    // cell
+    std::vector<double> cell(9, 0.0);
+    cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom;
+    cell[1] = ucell.latvec.e21 * ucell.lat0_angstrom;
+    cell[2] = ucell.latvec.e31 * ucell.lat0_angstrom;
+    cell[3] = ucell.latvec.e12 * ucell.lat0_angstrom;
+    cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom;
+    cell[5] = ucell.latvec.e32 * ucell.lat0_angstrom;
+    cell[6] = ucell.latvec.e13 * ucell.lat0_angstrom;
+    cell[7] = ucell.latvec.e23 * ucell.lat0_angstrom;
+    cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom;
+
+    // coord
+    std::vector<double> coord(3 * ucell.nat, 0.0);
+    int iat = 0;
+    const int nat = ucell.nat;
+    for (int it = 0; it < ucell.ntype; ++it)
+    {
+        for (int ia = 0; ia < ucell.atoms[it].na; ++ia)
+        {
+            coord[iat] = ucell.atoms[it].tau[ia].x * ucell.lat0_angstrom;
+            coord[iat + nat] = ucell.atoms[it].tau[ia].y * ucell.lat0_angstrom;
+            coord[iat + 2 * nat] = ucell.atoms[it].tau[ia].z * ucell.lat0_angstrom;
+            iat++;
+        }
+    }
+    assert(ucell.nat == iat);
+
+#ifdef __NEP
+    nep_potential = 0.0;
+    nep_force.zero_out();
+    nep_virial.zero_out();
+
+    nep.compute(atype, cell, coord, _e, _f, _v);
+
+    // unit conversion
+    const double fact_e = 1.0 / ModuleBase::Ry_to_eV;
+    const double fact_f = 1.0 / (ModuleBase::Ry_to_eV * ModuleBase::ANGSTROM_AU);
+    const double fact_v = 1.0 / (ucell.omega * ModuleBase::Ry_to_eV);
+
+
+    // potential energy
+    nep_potential = fact_e * std::accumulate(_e.begin(), _e.end(), 0.0) ;
+    GlobalV::ofs_running << " #TOTAL ENERGY# " << std::setprecision(11) << nep_potential * ModuleBase::Ry_to_eV << " eV"
+                         << std::endl;
+
+    // forces
+    for (int i = 0; i < nat; ++i)
+    {
+        nep_force(i, 0) = _f[i] * fact_f;
+        nep_force(i, 1) = _f[i + nat] * fact_f;
+        nep_force(i, 2) = _f[i + 2 * nat] * fact_f;
+    }
+
+    // virial
+    std::vector<double> v_sum(9, 0.0);
+    for (int j = 0; j < 9; ++j)
+    {
+        for (int i = 0; i < nat; ++i)
+        {
+            int index = j * nat + i;
+            v_sum[j] += _v[index];
+        }
+    }
+
+    // virial -> stress
+    for (int i = 0; i < 3; ++i)
+    {
+        for (int j = 0; j < 3; ++j)
+        {
+            nep_virial(i, j) = v_sum[3 * i + j] * fact_v;
+        }
+    }
+#else
+    ModuleBase::WARNING_QUIT("ESolver_NEP", "Please recompile with -D__NEP");
+#endif
+    ModuleBase::timer::tick("ESolver_NEP", "runner");
+}
+
+double ESolver_NEP::cal_energy()
+{
+    return nep_potential;
+}
+
+void ESolver_NEP::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
+{
+    force = nep_force;
+    ModuleIO::print_force(GlobalV::ofs_running, ucell, "TOTAL-FORCE (eV/Angstrom)", force, false);
+}
+
+void ESolver_NEP::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
+{
+    stress = nep_virial;
+    ModuleIO::print_stress("TOTAL-STRESS", stress, true, false, GlobalV::ofs_running);
+
+    // external stress
+    double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8;
+    double external_stress[3] = {PARAM.inp.press1, PARAM.inp.press2, PARAM.inp.press3};
+    for (int i = 0; i < 3; i++)
+    {
+        stress(i, i) -= external_stress[i] / unit_transform;
+    }
+}
+
+void ESolver_NEP::after_all_runners(UnitCell& ucell)
+{
+    GlobalV::ofs_running << "\n --------------------------------------------" << std::endl;
+    GlobalV::ofs_running << std::setprecision(16);
+    GlobalV::ofs_running << " !FINAL_ETOT_IS " << nep_potential * ModuleBase::Ry_to_eV << " eV" << std::endl;
+    GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
+}
+
+#ifdef __NEP
+void ESolver_NEP::type_map(const UnitCell& ucell)
+{   
+    // parse the element list from NEP model file
+    std::unordered_map<std::string, int> label;
+    std::string temp;
+    for (int i = 0; i < nep.element_list.size(); ++i)
+    {
+        label[nep.element_list[i]] = i; //> label: map from element string to index int.
+    }
+
+    std::cout << "\n Element list of model file " << nep_file << " " << std::endl;
+    std::cout << " ----------------------------------------------------------------";
+    int count = 0;
+    for (auto it = label.begin(); it != label.end(); ++it)
+    {
+        if (count % 5 == 0)
+        {
+            std::cout << std::endl;
+            std::cout << "  ";
+        }
+        count++;
+        temp = it->first + ": " + std::to_string(it->second);
+        std::cout << std::left << std::setw(10) << temp;
+    }
+    std::cout << "\n -----------------------------------------------------------------" << std::endl;
+
+    // parse the atype based on the element list
+    int iat = 0;
+    for (int it = 0; it < ucell.ntype; ++it)
+    {
+        for (int ia = 0; ia < ucell.atoms[it].na; ++ia)
+        {
+            if (label.find(ucell.atoms[it].label) == label.end())
+            {
+                ModuleBase::WARNING_QUIT("ESolver_NEP",
+                                         "The label " + ucell.atoms[it].label + " is not found in the type map.");
+            }
+            atype[iat] = label[ucell.atoms[it].label];
+            iat++;
+        }
+    }
+    assert(ucell.nat == iat);
+}
+#endif