feat: add CPU SNAP potential support

docs/advanced/md.md

@@ -96,4 +96,25 @@ And the filename of DP model is specified by keyword [pot_file](./input_files/in
 
 ## 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.
+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.
+
+## SNAP
+
+CPU SNAP calculations use `esolver_type snap`, with the coefficient file in
+`pot_file` and the LAMMPS-style parameter file in `snap_param_file`. The
+implementation supports linear and quadratic coefficients, chemical channels,
+outer cutoff switching, and `switchinnerflag` with `sinner` and `dinner`.
+
+For `chemflag 1`, the chemical-channel axes of each flattened coefficient vector
+follow the element order in the STRU file. Coefficient blocks are matched by
+element name, so their order within `pot_file` is not significant. Reordering
+the elements in STRU without reordering the chemical-channel coefficients can
+change the potential. The values in `sinner` and `dinner` also follow STRU
+element order. This corresponds to the first-appearance order of elements in
+LAMMPS `pair_coeff` mappings.
+
+The CPU implementation accepts but does not use the LAMMPS GPU execution
+options `chunksize` and `parallelthresh`. Unknown SNAP parameter keywords are
+rejected. Each element-pair cutoff must exceed `rmin0`; a contributing neighbor
+at or below `rmin0` causes an error. The radial mapping parameter `rfac0` must
+satisfy `0 < rfac0 <= 1` so that its angle does not cross pi inside a pair cutoff.

source/Makefile.Objects

@@ -33,6 +33,7 @@ VPATH=./src_global:\
 ./source_base/module_out:\
 ./source_md:\
 ./source_md/module_plumed:\
+./source_md/module_snap:\
 ./source_basis/module_pw:\
 ./source_base/module_fft:\
 ./source_esolver:\
@@ -344,6 +345,7 @@ OBJS_ESOLVER=esolver_factory.o\
     esolver_lj.o\
     esolver_dp.o\
     esolver_nep.o\
+    esolver_snap.o\
     esolver_of.o\
     esolver_of_tddft.o\
     esolver_of_interface.o\
@@ -561,6 +563,11 @@ OBJS_MD=fire.o\
     nhchain.o\
     plumed_interface.o\
     run_md.o\
+    module_snap/snap_calculator.o\
+    module_snap/snap_descriptor.o\
+    module_snap/snap_derivative.o\
+    module_snap/snap_model.o\
+    module_snap/snap_cg.o\
     verlet.o\
 
 OBJS_NEIGHBOR=sltk_atom.o\

source/source_esolver/CMakeLists.txt

@@ -11,6 +11,7 @@ list(
   esolver_lj.cpp
   esolver_dp.cpp
   esolver_nep.cpp
+  esolver_snap.cpp
   esolver_of.cpp
   esolver_of_tddft.cpp
   esolver_of_interface.cpp

source/source_esolver/esolver_factory.cpp

@@ -5,6 +5,7 @@
 #include "esolver_dfpt_pw.h"
 #include "esolver_sdft_pw.h"
 #include "esolver_ks_pw_tddft.h"
+#include "source_base/global_variable.h"
 #include "source_base/module_device/device.h"
 #include "source_hamilt/module_xc/general_exx_info.h"
 #include "source_io/module_parameter/parameter.h"
@@ -23,10 +24,12 @@
 #endif
 #include "esolver_dp.h"
 #include "esolver_nep.h"
+#include "esolver_snap.h"
 #include "esolver_lj.h"
 #include "esolver_of.h"
 #include "esolver_of_tddft.h"
 
+#include <sstream>
 #include <stdexcept>
 
 namespace ModuleESolver
@@ -113,12 +116,17 @@ std::string determine_type(const Input_para& inp)
     {
         esolver_type = "nep_pot";
     }
+    else if (inp.esolver_type == "snap")
+    {
+        esolver_type = "snap_pot";
+    }
     else if (esolver_type == "none")
     {
         ModuleBase::WARNING_QUIT("ESolver", "No such esolver_type combined with basis_type");
     }
 
-    GlobalV::ofs_running << "\n #ENERGY SOLVER# " << esolver_type << std::endl;
+    std::ostream& running_log = GlobalV::ofs_running;
+    running_log << "\n #ENERGY SOLVER# " << esolver_type << std::endl;
 
     auto device_info = inp.device;
 
@@ -129,8 +137,14 @@ std::string determine_type(const Input_para& inp)
             c = std::toupper(c);
         }
     }
-    base_device::information::output_device_info(std::cout, inp.device);
-    base_device::information::output_device_info(GlobalV::ofs_running, inp.device);
+    std::ostringstream device_summary;
+    base_device::information::output_device_info(device_summary, inp.device);
+    const std::string summary = device_summary.str();
+    if (!summary.empty())
+    {
+        std::cout << summary;
+        running_log << summary;
+    }
     /***auto end_time = std::chrono::high_resolution_clock::now();
     auto duration = std::chrono::duration_cast<std::chrono::duration<double>>(end_time - start_time);
     std::cout << "hipGetDeviceInfo took " << duration.count() << " seconds" << std::endl;***/
@@ -142,7 +156,7 @@ ESolver* init_esolver(const Input_para& inp)
 {
     // determine type of esolver based on INPUT information
     const std::string esolver_type = determine_type(inp);
-    const bool gamma_only = PARAM.globalv.gamma_only_local;
+    const bool gamma_only = inp.gamma_only;
 #ifdef __LCAO
     const LRRuntime lr_runtime{GlobalV::ofs_running, GlobalV::ofs_warning, GlobalV::MY_RANK};
 #endif
@@ -384,11 +398,16 @@ ESolver* init_esolver(const Input_para& inp)
     }
     else if (esolver_type == "dp_pot")
     {
-        return new ESolver_DP(PARAM.mdp.pot_file);
+        return new ESolver_DP(inp.mdp.pot_file);
     }
     else if (esolver_type == "nep_pot")
     {
-        return new ESolver_NEP(PARAM.mdp.pot_file);
+        return new ESolver_NEP(inp.mdp.pot_file);
+    }
+    else if (esolver_type == "snap_pot")
+    {
+        const int rank = GlobalV::MY_RANK;
+        return new ESolver_SNAP(inp.mdp.pot_file, inp.mdp.snap_param_file, GlobalV::ofs_running, rank);
     }
     throw std::invalid_argument("esolver_type = " + std::string(esolver_type) + ". Wrong in " + std::string(__FILE__)
                                 + " line " + std::to_string(__LINE__));