feat: add CPU SNAP potential support

docs/advanced/md.md

@@ -96,4 +96,24 @@ 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.

source/Makefile.Objects

@@ -32,6 +32,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:\
@@ -336,6 +337,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\
@@ -536,6 +538,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_descriptor_derivative.o\
+    module_snap/snap_model.o\
+    module_snap/snap_clebsch_gordan.o\
     verlet.o\
 
 OBJS_NEIGHBOR=sltk_atom.o\

source/source_esolver/CMakeLists.txt

@@ -9,6 +9,7 @@ list(APPEND objects
     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,6 +24,7 @@
 #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"
@@ -113,6 +115,10 @@ 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");
@@ -130,7 +136,6 @@ std::string determine_type(const Input_para& inp)
         }
     }
     base_device::information::output_device_info(std::cout, inp.device);
-    base_device::information::output_device_info(GlobalV::ofs_running, inp.device);
     /***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 +147,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 +389,15 @@ 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")
+    {
+        return new ESolver_SNAP(inp.mdp.pot_file, inp.mdp.snap_param_file, GlobalV::ofs_running);
     }
     throw std::invalid_argument("esolver_type = " + std::string(esolver_type) + ". Wrong in " + std::string(__FILE__)
                                 + " line " + std::to_string(__LINE__));