diff --git a/docs/advanced/md.md b/docs/advanced/md.md index 25c01b1c8a4..98d17de43ba 100644 --- a/docs/advanced/md.md +++ b/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. \ No newline at end of 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. diff --git a/source/Makefile.Objects b/source/Makefile.Objects index b2361a72e26..64a5b386506 100644 --- a/source/Makefile.Objects +++ b/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_derivative.o\ + module_snap/snap_model.o\ + module_snap/snap_cg.o\ verlet.o\ OBJS_NEIGHBOR=sltk_atom.o\ diff --git a/source/source_esolver/CMakeLists.txt b/source/source_esolver/CMakeLists.txt index a463b96f176..3ccfe94f3bd 100644 --- a/source/source_esolver/CMakeLists.txt +++ b/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 diff --git a/source/source_esolver/esolver_factory.cpp b/source/source_esolver/esolver_factory.cpp index de52d69a84c..24c6c51edec 100644 --- a/source/source_esolver/esolver_factory.cpp +++ b/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>(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,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__)); diff --git a/source/source_esolver/esolver_snap.cpp b/source/source_esolver/esolver_snap.cpp new file mode 100644 index 00000000000..a75321918ed --- /dev/null +++ b/source/source_esolver/esolver_snap.cpp @@ -0,0 +1,326 @@ +#include "esolver_snap.h" + +#include "source_base/constants.h" +#include "source_base/parallel_reduce.h" +#include "source_cell/mdcell.h" +#include "source_cell/module_neighlist/neighbor_search.h" +#include "source_io/module_output/output_log.h" +#include "source_io/module_parameter/parameter.h" + +#include +#include +#include +#include +#include + +namespace ModuleESolver +{ + +namespace +{ +void fail(const std::string& message, const int rank) +{ + if (rank != 0) + { + std::cerr << "ESolver_SNAP rank " << rank << ": " << message << std::endl; + } + ModuleBase::WARNING_QUIT("ESolver_SNAP", message); +} + +std::string full_path(const std::string& path, const int rank) +{ + char* resolved = ::realpath(path.c_str(), nullptr); + if (resolved == nullptr) + { + const std::string message = "Cannot resolve SNAP file path: " + path; + fail(message, rank); + } + const std::string result(resolved); + std::free(resolved); + return result; +} + +ModuleMD::SNAP::Neighbor make_neighbor(const NeighborAtom& center, + const NeighborAtom& neighbor, + const double lat0, + const ModuleMD::SNAP::ModelParameters& parameters) +{ + ModuleMD::SNAP::Neighbor result; + result.x = (neighbor.position_x - center.position_x) * lat0 * ModuleBase::BOHR_TO_A; + result.y = (neighbor.position_y - center.position_y) * lat0 * ModuleBase::BOHR_TO_A; + result.z = (neighbor.position_z - center.position_z) * lat0 * ModuleBase::BOHR_TO_A; + result.weight = parameters.element_weights[static_cast(neighbor.atom_type)]; + result.cutoff_radius = (parameters.element_radii[static_cast(center.atom_type)] + + parameters.element_radii[static_cast(neighbor.atom_type)]) + * parameters.rcutfac; + result.element_index = static_cast(neighbor.atom_type); + if (parameters.switchinner_flag) + { + const std::size_t center_type = static_cast(center.atom_type); + const std::size_t neighbor_type = static_cast(neighbor.atom_type); + result.sinner = 0.5 * (parameters.sinner[center_type] + parameters.sinner[neighbor_type]); + result.dinner = 0.5 * (parameters.dinner[center_type] + parameters.dinner[neighbor_type]); + } + return result; +} +} + +ESolver_SNAP::ESolver_SNAP(const std::string& coefficient_file, + const std::string& parameter_file, + std::ofstream& running_log, + const int rank) + : coefficient_file_(coefficient_file), parameter_file_(parameter_file), running_log_(running_log), rank_(rank) +{ + classname = "ESolver_SNAP"; +} + +void ESolver_SNAP::initialize(const std::vector& labels) +{ + std::string error; + if (!model_.load(coefficient_file_, parameter_file_, labels, error)) + { + fail(error, rank_); + } + const std::string coefficient_path = full_path(coefficient_file_, rank_); + const std::string parameter_path = full_path(parameter_file_, rank_); + running_log_ << " Read SNAP coefficient file: " << coefficient_path << std::endl; + running_log_ << " Read SNAP parameter file: " << parameter_path << std::endl; + if (rank_ == 0) + { + std::cout << " Read SNAP coefficient file: " << coefficient_path << std::endl; + std::cout << " Read SNAP parameter file: " << parameter_path << std::endl; + } + try + { + calculator_.reset(new ModuleMD::SNAP::Calculator(model_)); + } + catch (const std::exception& error) + { + const std::string message = std::string("Invalid SNAP model: ") + error.what(); + fail(message, rank_); + } + snap_virial_.create(3, 3); +} + +void ESolver_SNAP::before_all_runners(BaseCell& cell, const Input_para& inp) +{ + inp_ = &inp; + if (cell.kind() == BaseCell::Kind::mdcell) + { + MDCell& mdcell = static_cast(cell); + initialize(mdcell.type_labels()); + mdcell.set_neighbor_cutoff(model_.parameters.rcut * ModuleBase::ANGSTROM_AU); + return; + } + + UnitCell& ucell = static_cast(cell); + std::vector labels; + for (int it = 0; it < ucell.ntype; ++it) + { + labels.push_back(ucell.atoms[it].label); + } + initialize(labels); + snap_force_.create(ucell.nat, 3); +} + +void ESolver_SNAP::runner(BaseCell& cell, const int istep) try +{ + static_cast(istep); + if (!calculator_) + { + fail("SNAP solver was not initialized.", rank_); + } + snap_potential_ = 0.0; + snap_virial_.zero_out(); + + if (cell.kind() == BaseCell::Kind::mdcell) + { + MDCell& mdcell = static_cast(cell); + if (!mdcell.has_neighbor_search()) + { + fail("MDCell neighbors must be prepared before ESolver_SNAP::runner().", rank_); + } + std::vector& owned_atoms = mdcell.owned_atoms(); + const std::vector& all_atoms = mdcell.neighbor_search().get_all_atoms(); + const NeighborList& list = mdcell.neighbor_search().get_neighbor_list(); + const int nowned = static_cast(owned_atoms.size()); + ModuleMD::SNAP::CellAccumulator contribution(all_atoms.size()); + for (int i = 0; i < nowned; ++i) + { + const NeighborAtom center(owned_atoms[static_cast(i)].cart.x, + owned_atoms[static_cast(i)].cart.y, + owned_atoms[static_cast(i)].cart.z, + owned_atoms[static_cast(i)].type, + i, + i); + std::vector neighbors; + std::vector neighbor_indices; + for (int k = 0; k < list.get_numneigh(i); ++k) + { + const int neighbor_index = list.get_firstneigh(i)[k]; + const NeighborAtom& neighbor = all_atoms[static_cast(neighbor_index)]; + if (neighbor.atom_type < 0 || static_cast(neighbor.atom_type) >= model_.parameters.elements.size()) + { + fail("MDCell atom type is outside the SNAP element map.", rank_); + } + neighbors.push_back(make_neighbor(center, neighbor, mdcell.lat0(), model_.parameters)); + neighbor_indices.push_back(static_cast(neighbor_index)); + } + const std::size_t center_element = static_cast(center.atom_type); + const ModuleMD::SNAP::ForceResult result = calculator_->compute_with_analytic_force(neighbors, center_element); + const std::size_t center_index = static_cast(i); + contribution.add(center_index, neighbor_indices, neighbors, result); + } + snap_potential_ = contribution.energy() / ModuleBase::Ry_to_eV; + const double virial_scale = 1.0 / (ModuleBase::Ry_to_eV * mdcell.omega()); + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + const std::size_t index = static_cast(3 * row + column); + snap_virial_(row, column) = contribution.virial()[index] * virial_scale; + } + } + const double force_scale = 1.0 / (ModuleBase::Ry_to_eV * ModuleBase::ANGSTROM_AU); + for (std::size_t i = 0; i < mdcell.ghost_atoms().size(); ++i) + { + const std::array& value = contribution.forces()[nowned + i]; + const double fx = value[0] * force_scale; + const double fy = value[1] * force_scale; + const double fz = value[2] * force_scale; + mdcell.ghost_atoms()[i].force.set(fx, fy, fz); + } + Parallel_Reduce::reduce_all(snap_potential_); + Parallel_Reduce::reduce_all(snap_virial_.c, 9); + for (int i = 0; i < nowned; ++i) + { + const std::array& value = contribution.forces()[static_cast(i)]; + const double fx = value[0] * force_scale; + const double fy = value[1] * force_scale; + const double fz = value[2] * force_scale; + owned_atoms[static_cast(i)].force.set(fx, fy, fz); + } + return; + } + + UnitCell& ucell = static_cast(cell); + NeighborSearch search; + search.init(ucell, model_.parameters.rcut * ModuleBase::ANGSTROM_AU); + search.build_neighbors(); + const NeighborList& list = search.get_neighbor_list(); + const std::vector& inside = search.get_inside_atoms(); + const std::vector& all = search.get_all_atoms(); + snap_force_.zero_out(); + const std::size_t atom_count = static_cast(ucell.nat); + ModuleMD::SNAP::CellAccumulator contribution(atom_count); + std::vector offsets(static_cast(ucell.ntype + 1), 0); + for (int it = 0; it < ucell.ntype; ++it) + { + offsets[static_cast(it + 1)] = offsets[static_cast(it)] + ucell.atoms[it].na; + } + for (int local_i = 0; local_i < list.get_ncentral_atoms(); ++local_i) + { + const NeighborAtom& center = inside[static_cast(local_i)]; + std::vector neighbors; + std::vector atom_indices; + for (int k = 0; k < list.get_numneigh(local_i); ++k) + { + const NeighborAtom& neighbor = all[static_cast(list.get_firstneigh(local_i)[k])]; + neighbors.push_back(make_neighbor(center, neighbor, ucell.lat0, model_.parameters)); + const int atom_index = offsets[static_cast(neighbor.atom_type)] + neighbor.atom_index; + const std::size_t index = static_cast(atom_index); + atom_indices.push_back(index); + } + const std::size_t center_element = static_cast(center.atom_type); + const ModuleMD::SNAP::ForceResult result = calculator_->compute_with_analytic_force(neighbors, center_element); + const int center_index = offsets[static_cast(center.atom_type)] + center.atom_index; + const std::size_t index = static_cast(center_index); + contribution.add(index, atom_indices, neighbors, result); + } + snap_potential_ = contribution.energy() / ModuleBase::Ry_to_eV; + const double virial_scale = 1.0 / (ModuleBase::Ry_to_eV * ucell.omega); + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + const std::size_t index = static_cast(3 * row + column); + snap_virial_(row, column) = contribution.virial()[index] * virial_scale; + } + } + const double force_scale = 1.0 / (ModuleBase::Ry_to_eV * ModuleBase::ANGSTROM_AU); + for (int atom_index = 0; atom_index < ucell.nat; ++atom_index) + { + for (int direction = 0; direction < 3; ++direction) + { + const std::size_t index = static_cast(atom_index); + const std::size_t component = static_cast(direction); + snap_force_(atom_index, direction) = contribution.forces()[index][component] * force_scale; + } + } +} +catch (const std::exception& error) +{ + const std::string message = std::string("SNAP evaluation failed: ") + error.what(); + fail(message, rank_); +} + +double ESolver_SNAP::cal_energy() +{ + return snap_potential_; +} + +void ESolver_SNAP::cal_force(BaseCell& cell, ModuleBase::matrix& force) +{ + if (cell.kind() == BaseCell::Kind::mdcell) + { + const MDCell& mdcell = static_cast(cell); + force.create(mdcell.owned_atoms().size(), 3); + for (std::size_t i = 0; i < mdcell.owned_atoms().size(); ++i) + { + force(i, 0) = mdcell.owned_atoms()[i].force.x; + force(i, 1) = mdcell.owned_atoms()[i].force.y; + force(i, 2) = mdcell.owned_atoms()[i].force.z; + } + return; + } + force = snap_force_; + ModuleIO::print_force(running_log_, + static_cast(cell), + "TOTAL-FORCE (eV/Angstrom)", + force, + false); +} + +void ESolver_SNAP::cal_stress(BaseCell& cell, ModuleBase::matrix& stress) +{ + stress = snap_virial_; + + if (cell.kind() == BaseCell::Kind::unitcell) + { + ModuleIO::print_stress("TOTAL-STRESS", stress, true, false, running_log_); + } + else + { +#ifdef __MPI + if (static_cast(cell).mpi_rank() == 0) +#endif + { + ModuleIO::print_stress("TOTAL-STRESS", stress, true, false, running_log_); + } + } + + const double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; + const double external_stress[3] = {inp_->press1, inp_->press2, inp_->press3}; + for (int i = 0; i < 3; ++i) + { + stress(i, i) -= external_stress[i] / unit_transform; + } +} + +void ESolver_SNAP::after_all_runners(BaseCell& cell) +{ + static_cast(cell); + running_log_ << "\n !FINAL_ETOT_IS " << snap_potential_ * ModuleBase::Ry_to_eV << " eV\n"; +} + +} // namespace ModuleESolver diff --git a/source/source_esolver/esolver_snap.h b/source/source_esolver/esolver_snap.h new file mode 100644 index 00000000000..6206750b201 --- /dev/null +++ b/source/source_esolver/esolver_snap.h @@ -0,0 +1,46 @@ +#ifndef ESOLVER_SNAP_H +#define ESOLVER_SNAP_H + +#include "esolver.h" +#include "source_md/module_snap/snap_calculator.h" + +#include +#include +#include +#include + +namespace ModuleESolver +{ + +class ESolver_SNAP : public ESolver +{ + public: + ESolver_SNAP(const std::string& coefficient_file, + const std::string& parameter_file, + std::ofstream& running_log, + int rank); + + void before_all_runners(BaseCell& cell, const Input_para& inp) override; + void runner(BaseCell& cell, int istep) override; + double cal_energy() override; + void cal_force(BaseCell& cell, ModuleBase::matrix& force) override; + void cal_stress(BaseCell& cell, ModuleBase::matrix& stress) override; + void after_all_runners(BaseCell& cell) override; + + private: + void initialize(const std::vector& labels); + + std::string coefficient_file_; + std::string parameter_file_; + std::ofstream& running_log_; + const int rank_; + ModuleMD::SNAP::ModelData model_; + std::unique_ptr calculator_; + double snap_potential_ = 0.0; + ModuleBase::matrix snap_force_; + ModuleBase::matrix snap_virial_; +}; + +} // namespace ModuleESolver + +#endif // ESOLVER_SNAP_H diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 75be1fba4b8..ff9b3771183 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -21,7 +21,7 @@ struct Input_para ///< "nscf" : non-self consistent calculation. ///< "relax" : cell relaxations bool socket_driver = false; ///< run ABACUS as an i-PI socket client - std::string esolver_type = "ksdft"; ///< the energy solver: ksdft, sdft, ofdft, tddft, lj, dp + std::string esolver_type = "ksdft"; ///< the energy solver: ksdft, sdft, ofdft, tddft, lj, dp, nep, snap /* symmetry level: -1, no symmetry at all; 0, only basic time reversal would be considered; diff --git a/source/source_io/module_parameter/md_parameter.h b/source/source_io/module_parameter/md_parameter.h index 706218e6960..cdc601c8502 100644 --- a/source/source_io/module_parameter/md_parameter.h +++ b/source/source_io/module_parameter/md_parameter.h @@ -31,6 +31,7 @@ struct MD_para std::vector lj_epsilon = {}; ///< the value of epsilon for LJ potential (eV) std::vector lj_sigma = {}; ///< the value of sigma for LJ potential (\AA) std::string pot_file = "graph.pb"; ///< the filename of potential files for CMD such as DP + std::string snap_param_file = "snap.param"; ///< SNAP parameter file double dp_rescaling = 1.0; ///< rescaling factor for DP model. Energy, force and stress will be multiplied by this factor. std::vector dp_fparam = {}; ///< The frame parameter for dp potential. The array can be of size : diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 43adbfc8169..535bde992cf 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -247,6 +247,7 @@ Socket mode always computes energy. Force and stress extraction follows cal_forc * `lj`: Lennard-Jones potential. * `dp`: Deep Potential (DP) model. * `nep`: Neuroevolution potential (NEP). +* `snap`: Spectral Neighbor Analysis Potential (SNAP). * `ks-lr`: Kohn-Sham density functional theory followed by linear-response time-dependent density functional theory (LR-TDDFT; under development). * `lr`: Linear-response time-dependent density functional theory (LR-TDDFT) using precomputed Kohn-Sham orbitals (under development). * `dfpt`: Density functional perturbation theory (DFPT; under development). @@ -255,7 +256,7 @@ Socket mode always computes energy. Force and stress extraction follows cal_forc item.default_value = "ksdft"; read_sync_string(input.esolver_type); item.check_value = [](const Input_Item& item, const Parameter& para) { - const std::vector esolver_types = { "ksdft", "sdft", "ofdft", "tdofdft", "tddft", "lj", "dp", "nep", "lr", "ks-lr", "dfpt" }; + const std::vector esolver_types = { "ksdft", "sdft", "ofdft", "tdofdft", "tddft", "lj", "dp", "nep", "snap", "lr", "ks-lr", "dfpt" }; if (std::find(esolver_types.begin(), esolver_types.end(), para.input.esolver_type) == esolver_types.end()) { const std::string warningstr = nofound_str(esolver_types, "esolver_type"); @@ -265,12 +266,17 @@ Socket mode always computes energy. Force and stress extraction follows cal_forc { ModuleBase::WARNING_QUIT("ReadInput", "esolver_type=tddft requires calculation=md."); } - if (para.input.esolver_type == "dp" || para.input.esolver_type == "nep") + if (para.input.esolver_type == "dp" || para.input.esolver_type == "nep" + || para.input.esolver_type == "snap") { if (access(para.input.mdp.pot_file.c_str(), 0) == -1) { ModuleBase::WARNING_QUIT("ReadInput", "Can not find `pot_file` !"); } + if (para.input.esolver_type == "snap" && access(para.input.mdp.snap_param_file.c_str(), 0) == -1) + { + ModuleBase::WARNING_QUIT("ReadInput", "Can not find `snap_param_file` !"); + } } // LR reads the ground state wave function from a separate SCF run, // so it cannot be combined with a self-consistent calculation. 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 a41343fd261..2602e1c421e 100644 --- a/source/source_io/module_parameter/read_input_item_md.cpp +++ b/source/source_io/module_parameter/read_input_item_md.cpp @@ -552,6 +552,18 @@ Note: It is a system-dependent empirical parameter. An improper choice might lea read_sync_string(input.mdp.pot_file); this->add_item(item); } + { + Input_Item item("snap_param_file"); + item.annotation = "the SNAP parameter file"; + item.category = "Molecular dynamics"; + item.type = "String"; + item.description = "The SNAP parameter file used with pot_file as the SNAP coefficient file."; + item.default_value = "snap.param"; + item.unit = ""; + item.set_availability("esolver_type==snap"); + read_sync_string(input.mdp.snap_param_file); + this->add_item(item); + } { Input_Item item("dp_rescaling"); item.annotation = "rescaling factor for dp potential"; diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 67da741a041..e41720119bd 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -52,7 +52,8 @@ void Driver::driver_run() this->init_hardware(); ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp); - const bool direct_mdcell = input.esolver_type == "lj" || input.esolver_type == "dp" || input.esolver_type == "nep"; + const bool direct_mdcell = input.esolver_type == "lj" || input.esolver_type == "dp" || input.esolver_type == "nep" + || input.esolver_type == "snap"; UnitCell ucell; bool ucell_initialized = false; const auto initialize_ucell = [&ucell, &ucell_initialized, &input]() { diff --git a/source/source_md/CMakeLists.txt b/source/source_md/CMakeLists.txt index ac5d35156ec..b71246d8370 100644 --- a/source/source_md/CMakeLists.txt +++ b/source/source_md/CMakeLists.txt @@ -9,6 +9,8 @@ add_library( run_md.cpp verlet.cpp) +add_subdirectory(module_snap) + if(ENABLE_COVERAGE) add_coverage(md) endif() diff --git a/source/source_md/module_snap/CMakeLists.txt b/source/source_md/module_snap/CMakeLists.txt new file mode 100644 index 00000000000..c10bf8ed1d0 --- /dev/null +++ b/source/source_md/module_snap/CMakeLists.txt @@ -0,0 +1,13 @@ +target_sources( + md + PRIVATE + snap_calculator.cpp + snap_descriptor.cpp + snap_derivative.cpp + snap_model.cpp + snap_cg.cpp +) + +if(BUILD_TESTING) + add_subdirectory(unittest) +endif() diff --git a/source/source_md/module_snap/snap_calculator.cpp b/source/source_md/module_snap/snap_calculator.cpp new file mode 100644 index 00000000000..86440767fd7 --- /dev/null +++ b/source/source_md/module_snap/snap_calculator.cpp @@ -0,0 +1,273 @@ +#include "snap_calculator.h" + +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +Calculator::Calculator(const ModelData& model) + : model_(model), + descriptor_evaluator_(model.parameters), + descriptor_size_(0) +{ + const std::size_t channels = descriptor_evaluator_.channel_count(); + descriptor_size_ = descriptor_evaluator_.bispectrum_size() * channels * channels * channels; + if (model_.coefficients.size() != model_.parameters.elements.size()) + { + throw std::invalid_argument("SNAP coefficient blocks do not match the element list."); + } + const std::size_t expected_count = descriptor_size_ + 1 + + (model_.parameters.quadratic_flag + ? descriptor_size_ * (descriptor_size_ + 1) / 2 + : 0); + for (std::size_t element = 0; element < model_.coefficients.size(); ++element) + { + if (model_.coefficients[element].size() != expected_count) + { + throw std::invalid_argument("SNAP coefficient count for element " + model_.parameters.elements[element] + + " does not match the descriptor size: expected " + + std::to_string(expected_count) + ", found " + + std::to_string(model_.coefficients[element].size()) + "."); + } + } + if (model_.parameters.bzero_flag) + { + const std::vector empty_neighbors; + bzero_.reserve(model_.coefficients.size()); + for (std::size_t element = 0; element < model_.coefficients.size(); ++element) + { + const std::vector > self_u + = descriptor_evaluator_.accumulate_u(empty_neighbors, element); + const std::vector self_bispectrum = descriptor_evaluator_.compute_bispectrum(self_u); + bzero_.push_back(self_bispectrum); + } + } +} + +EnergyResult Calculator::compute(const std::vector& neighbors, const std::size_t central_element) const +{ + if (central_element >= model_.coefficients.size()) + { + throw std::out_of_range("SNAP central element index is outside the model."); + } + + const std::vector > u = descriptor_evaluator_.accumulate_u(neighbors, central_element); + return compute_from_u(u, central_element); +} + +EnergyResult Calculator::compute_from_u(const std::vector >& u, + const std::size_t central_element) const +{ + EnergyResult result; + result.descriptor = descriptor_evaluator_.compute_bispectrum(u); + if (model_.parameters.bzero_flag) + { + for (std::size_t index = 0; index < result.descriptor.size(); ++index) + { + result.descriptor[index] -= bzero_[central_element][index]; + } + } + const std::vector& coefficients = model_.coefficients[central_element]; + result.energy = coefficients[0]; + for (std::size_t index = 0; index < result.descriptor.size(); ++index) + { + result.energy += coefficients[index + 1] * result.descriptor[index]; + } + if (model_.parameters.quadratic_flag) + { + std::size_t coefficient = result.descriptor.size() + 1; + for (std::size_t i = 0; i < result.descriptor.size(); ++i) + { + for (std::size_t j = i; j < result.descriptor.size(); ++j) + { + result.energy += coefficients[coefficient++] * result.descriptor[i] * result.descriptor[j] + * (i == j ? 0.5 : 1.0); + } + } + } + return result; +} + +ForceResult Calculator::compute_with_analytic_force(const std::vector& neighbors, + const std::size_t central_element) const +{ + if (central_element >= model_.coefficients.size()) + { + throw std::out_of_range("SNAP central element index is outside the model."); + } + ForceResult result; + const std::vector > u = descriptor_evaluator_.accumulate_u(neighbors, central_element); + result.energy = compute_from_u(u, central_element); + result.center_force.assign(3, 0.0); + result.neighbor_forces.assign(neighbors.size(), std::vector(3, 0.0)); + const std::vector& coefficients = model_.coefficients[central_element]; + std::vector force_coefficients(result.energy.descriptor.size(), 0.0); + for (std::size_t i = 0; i < force_coefficients.size(); ++i) + { + force_coefficients[i] = coefficients[i + 1]; + } + if (model_.parameters.quadratic_flag) + { + std::size_t coefficient = force_coefficients.size() + 1; + for (std::size_t i = 0; i < force_coefficients.size(); ++i) + { + for (std::size_t j = i; j < force_coefficients.size(); ++j) + { + const double value = coefficients[coefficient++]; + if (i == j) + { + force_coefficients[i] += value * result.energy.descriptor[i]; + } + else + { + force_coefficients[i] += value * result.energy.descriptor[j]; + force_coefficients[j] += value * result.energy.descriptor[i]; + } + } + } + } + for (std::size_t neighbor_index = 0; neighbor_index < neighbors.size(); ++neighbor_index) + { + for (int direction = 0; direction < 3; ++direction) + { + const UWithDerivative u_derivative = + descriptor_evaluator_.compute_u_with_derivative(neighbors[neighbor_index], direction); + std::vector > du(u.size(), std::complex(0.0, 0.0)); + const std::size_t derivative_offset = model_.parameters.chem_flag + ? neighbors[neighbor_index].element_index + * descriptor_evaluator_.u_size() + : 0; + const std::vector >& one_neighbor_derivative = u_derivative.derivatives; + for (std::size_t index = 0; index < one_neighbor_derivative.size(); ++index) + { + du[derivative_offset + index] = one_neighbor_derivative[index]; + } + const std::vector db = descriptor_evaluator_.compute_bispectrum_derivative(u, du); + double derivative = 0.0; + for (std::size_t index = 0; index < db.size(); ++index) + { + derivative += force_coefficients[index] * db[index]; + } + result.center_force[static_cast(direction)] += derivative; + result.neighbor_forces[neighbor_index][static_cast(direction)] = -derivative; + } + } + return result; +} + +ForceResult Calculator::compute_with_finite_difference_force(const std::vector& neighbors, + const std::size_t central_element, + const double displacement) const +{ + if (displacement <= 0.0) + { + throw std::invalid_argument("SNAP finite-difference displacement must be positive."); + } + ForceResult result; + result.energy = compute(neighbors, central_element); + result.center_force.assign(3, 0.0); + result.neighbor_forces.assign(neighbors.size(), std::vector(3, 0.0)); + + for (std::size_t neighbor_index = 0; neighbor_index < neighbors.size(); ++neighbor_index) + { + for (int direction = 0; direction < 3; ++direction) + { + std::vector plus = neighbors; + std::vector minus = neighbors; + if (direction == 0) + { + plus[neighbor_index].x += displacement; + minus[neighbor_index].x -= displacement; + } + else if (direction == 1) + { + plus[neighbor_index].y += displacement; + minus[neighbor_index].y -= displacement; + } + else + { + plus[neighbor_index].z += displacement; + minus[neighbor_index].z -= displacement; + } + const double derivative = + (compute(plus, central_element).energy - compute(minus, central_element).energy) / (2.0 * displacement); + result.center_force[static_cast(direction)] += derivative; + result.neighbor_forces[neighbor_index][static_cast(direction)] -= derivative; + } + } + return result; +} + +CellAccumulator::CellAccumulator(const std::size_t atom_count) : forces_(atom_count) +{ + for (std::size_t atom = 0; atom < forces_.size(); ++atom) + { + forces_[atom].fill(0.0); + } + virial_.fill(0.0); +} + +void CellAccumulator::add(const std::size_t center_index, + const std::vector& neighbor_indices, + const std::vector& neighbors, + const ForceResult& result) +{ + if (center_index >= forces_.size() || neighbor_indices.size() != neighbors.size() + || result.neighbor_forces.size() != neighbors.size() || result.center_force.size() != 3) + { + throw std::invalid_argument("SNAP cell contribution has inconsistent atom or force indices."); + } + for (std::size_t index = 0; index < neighbor_indices.size(); ++index) + { + if (neighbor_indices[index] >= forces_.size() || result.neighbor_forces[index].size() != 3) + { + throw std::invalid_argument("SNAP neighbor force index is outside the cell."); + } + } + + energy_ += result.energy.energy; + for (int direction = 0; direction < 3; ++direction) + { + forces_[center_index][static_cast(direction)] += + result.center_force[static_cast(direction)]; + } + for (std::size_t index = 0; index < neighbors.size(); ++index) + { + const Neighbor& neighbor = neighbors[index]; + const double displacement[3] = {neighbor.x, neighbor.y, neighbor.z}; + const std::size_t atom = neighbor_indices[index]; + for (int direction = 0; direction < 3; ++direction) + { + forces_[atom][static_cast(direction)] += + result.neighbor_forces[index][static_cast(direction)]; + } + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + virial_[static_cast(3 * row + column)] += + displacement[row] * result.neighbor_forces[index][static_cast(column)]; + } + } + } +} + +double CellAccumulator::energy() const +{ + return energy_; +} + +const std::vector >& CellAccumulator::forces() const +{ + return forces_; +} + +const std::array& CellAccumulator::virial() const +{ + return virial_; +} + +} // namespace SNAP +} // namespace ModuleMD diff --git a/source/source_md/module_snap/snap_calculator.h b/source/source_md/module_snap/snap_calculator.h new file mode 100644 index 00000000000..077cd4f5ca3 --- /dev/null +++ b/source/source_md/module_snap/snap_calculator.h @@ -0,0 +1,75 @@ +#ifndef SNAP_CALCULATOR_H +#define SNAP_CALCULATOR_H + +#include "snap_descriptor.h" +#include "snap_model.h" + +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +struct EnergyResult +{ + double energy = 0.0; + std::vector descriptor; +}; + +struct ForceResult +{ + EnergyResult energy; + std::vector center_force; + std::vector > neighbor_forces; +}; + +class Calculator +{ + public: + explicit Calculator(const ModelData& model); + + EnergyResult compute(const std::vector& neighbors, std::size_t central_element) const; + + ForceResult compute_with_analytic_force(const std::vector& neighbors, + std::size_t central_element) const; + + ForceResult compute_with_finite_difference_force(const std::vector& neighbors, + std::size_t central_element, + double displacement) const; + + private: + EnergyResult compute_from_u(const std::vector >& u, + std::size_t central_element) const; + + ModelData model_; + DescriptorEvaluator descriptor_evaluator_; + std::size_t descriptor_size_; + std::vector > bzero_; +}; + +class CellAccumulator +{ + public: + explicit CellAccumulator(std::size_t atom_count); + + void add(std::size_t center_index, + const std::vector& neighbor_indices, + const std::vector& neighbors, + const ForceResult& result); + + double energy() const; + const std::vector >& forces() const; + const std::array& virial() const; + + private: + double energy_ = 0.0; + std::vector > forces_; + std::array virial_; +}; + +} // namespace SNAP +} // namespace ModuleMD + +#endif // SNAP_CALCULATOR_H diff --git a/source/source_md/module_snap/snap_cg.cpp b/source/source_md/module_snap/snap_cg.cpp new file mode 100644 index 00000000000..e9f2ca0ffb7 --- /dev/null +++ b/source/source_md/module_snap/snap_cg.cpp @@ -0,0 +1,153 @@ +#include "snap_cg.h" + +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +namespace +{ +int half_sum(const int left, const int right) +{ + return (left + right) / 2; +} +} + +ClebschGordanTable::ClebschGordanTable(const int twojmax) +{ + twojmax_ = twojmax < 0 ? 0 : twojmax; + for (int j1 = 0; j1 <= twojmax_; ++j1) + { + for (int j2 = 0; j2 <= twojmax_; ++j2) + { + for (int j = std::abs(j1 - j2); j <= std::min(twojmax_, j1 + j2); j += 2) + { + for (int m1 = -j1; m1 <= j1; m1 += 2) + { + for (int m2 = -j2; m2 <= j2; m2 += 2) + { + const int m = m1 + m2; + if (m < -j || m > j) + { + continue; + } + Entry entry; + entry.j1 = j1; + entry.m1 = m1; + entry.j2 = j2; + entry.m2 = m2; + entry.j = j; + entry.m = m; + entry.value = coefficient_uncached(j1, m1, j2, m2, j, m); + const std::tuple key(j1, m1, j2, m2, j, m); + entry_indices_[key] = entries_.size(); + entries_.push_back(entry); + } + } + } + } + } +} + +double ClebschGordanTable::factorial(const int value) +{ + if (value < 0) + { + return 0.0; + } + return std::tgamma(static_cast(value) + 1.0); +} + +bool ClebschGordanTable::allowed(const int j1, + const int m1, + const int j2, + const int m2, + const int j, + const int m) +{ + if (j1 < 0 || j2 < 0 || j < 0 || m1 < -j1 || m1 > j1 || m2 < -j2 || m2 > j2 || m < -j || m > j) + { + return false; + } + if ((j1 - m1) % 2 != 0 || (j2 - m2) % 2 != 0 || (j - m) % 2 != 0) + { + return false; + } + if ((j1 + j2 + j) % 2 != 0 || j < std::abs(j1 - j2) || j > j1 + j2) + { + return false; + } + return m == m1 + m2; +} + +double ClebschGordanTable::coefficient_uncached(const int j1, + const int m1, + const int j2, + const int m2, + const int j, + const int m) +{ + if (!allowed(j1, m1, j2, m2, j, m)) + { + return 0.0; + } + + const double delta = std::sqrt( + factorial((j1 + j2 - j) / 2) * factorial((j1 - j2 + j) / 2) + * factorial((-j1 + j2 + j) / 2) / factorial((j1 + j2 + j) / 2 + 1)); + const double prefactor = std::sqrt( + (j + 1.0) * factorial(half_sum(j1, m1)) * factorial(half_sum(j1, -m1)) + * factorial(half_sum(j2, m2)) * factorial(half_sum(j2, -m2)) + * factorial(half_sum(j, m)) * factorial(half_sum(j, -m))); + + double sum = 0.0; + const int zmin = std::max(0, std::max(-(j - j2 + m1) / 2, -(j - j1 - m2) / 2)); + const int zmax = std::min((j1 + j2 - j) / 2, + std::min((j1 - m1) / 2, (j2 + m2) / 2)); + for (int z = zmin; z <= zmax; ++z) + { + const int a = z; + const int b = (j1 + j2 - j) / 2 - z; + const int c = (j1 - m1) / 2 - z; + const int d = (j2 + m2) / 2 - z; + const int e = (j - j2 + m1) / 2 + z; + const int f = (j - j1 - m2) / 2 + z; + if (a < 0 || b < 0 || c < 0 || d < 0 || e < 0 || f < 0) + { + continue; + } + const double denominator = factorial(a) * factorial(b) * factorial(c) * factorial(d) * factorial(e) * factorial(f); + sum += (z % 2 == 0 ? 1.0 : -1.0) / denominator; + } + + return delta * prefactor * sum; +} + +const ClebschGordanTable::Entry* ClebschGordanTable::find(const int j1, + const int m1, + const int j2, + const int m2, + const int j, + const int m) const +{ + const std::tuple key(j1, m1, j2, m2, j, m); + const std::map, std::size_t>::const_iterator found + = entry_indices_.find(key); + return found == entry_indices_.end() ? nullptr : &entries_[found->second]; +} + +double ClebschGordanTable::coefficient(const int j1, + const int m1, + const int j2, + const int m2, + const int j, + const int m) const +{ + const Entry* entry = find(j1, m1, j2, m2, j, m); + return entry == nullptr ? 0.0 : entry->value; +} + +} // namespace SNAP +} // namespace ModuleMD diff --git a/source/source_md/module_snap/snap_cg.h b/source/source_md/module_snap/snap_cg.h new file mode 100644 index 00000000000..7198536882c --- /dev/null +++ b/source/source_md/module_snap/snap_cg.h @@ -0,0 +1,51 @@ +#ifndef SNAP_CG_H +#define SNAP_CG_H + +#include +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +/** + * @brief Clebsch-Gordan coefficients indexed with doubled angular momenta. + * + * If the stored index is j, the physical angular momentum is j/2. The same + * convention is used for m and m'. + */ +class ClebschGordanTable +{ + public: + explicit ClebschGordanTable(int twojmax); + + double coefficient(int j1, int m1, int j2, int m2, int j, int m) const; + + private: + struct Entry + { + int j1; + int m1; + int j2; + int m2; + int j; + int m; + double value; + }; + + static bool allowed(int j1, int m1, int j2, int m2, int j, int m); + static double coefficient_uncached(int j1, int m1, int j2, int m2, int j, int m); + static double factorial(int value); + const Entry* find(int j1, int m1, int j2, int m2, int j, int m) const; + + int twojmax_; + std::vector entries_; + std::map, std::size_t> entry_indices_; +}; + +} // namespace SNAP +} // namespace ModuleMD + +#endif // SNAP_CG_H diff --git a/source/source_md/module_snap/snap_derivative.cpp b/source/source_md/module_snap/snap_derivative.cpp new file mode 100644 index 00000000000..0448560f056 --- /dev/null +++ b/source/source_md/module_snap/snap_derivative.cpp @@ -0,0 +1,190 @@ +#include "snap_descriptor.h" + +#include "source_base/constants.h" + +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +namespace +{ +struct Dual +{ + double value; + double derivative; + + Dual(const double value_in = 0.0, const double derivative_in = 0.0) + : value(value_in), derivative(derivative_in) + { + } +}; + +Dual operator+(const Dual& left, const Dual& right) +{ + return Dual(left.value + right.value, left.derivative + right.derivative); +} + +Dual operator-(const Dual& left, const Dual& right) +{ + return Dual(left.value - right.value, left.derivative - right.derivative); +} + +Dual operator-(const Dual& value) +{ + return Dual(-value.value, -value.derivative); +} + +Dual operator*(const Dual& left, const Dual& right) +{ + return Dual(left.value * right.value, + left.derivative * right.value + left.value * right.derivative); +} + +Dual operator/(const Dual& left, const Dual& right) +{ + const double inverse = 1.0 / right.value; + return Dual(left.value * inverse, + (left.derivative * right.value - left.value * right.derivative) * inverse * inverse); +} + +Dual operator+(const Dual& left, const double right) +{ + return left + Dual(right); +} + +Dual operator-(const Dual& left, const double right) +{ + return left - Dual(right); +} + +Dual operator*(const Dual& left, const double right) +{ + return left * Dual(right); +} + +Dual operator*(const double left, const Dual& right) +{ + return Dual(left) * right; +} + +Dual operator/(const Dual& left, const double right) +{ + return left / Dual(right); +} + +Dual operator/(const double left, const Dual& right) +{ + return Dual(left) / right; +} + +Dual& operator*=(Dual& left, const Dual& right) +{ + left = left * right; + return left; +} + +Dual& operator+=(Dual& left, const Dual& right) +{ + left = left + right; + return left; +} + +Dual dual_sqrt(const Dual& value) +{ + const double root = std::sqrt(value.value); + return Dual(root, value.derivative / (2.0 * root)); +} + +Dual dual_tan(const Dual& value) +{ + const double cosine = std::cos(value.value); + return Dual(std::tan(value.value), value.derivative / (cosine * cosine)); +} +} + +UWithDerivative DescriptorEvaluator::compute_u_with_derivative(const Neighbor& neighbor, const int direction) const +{ + if (direction < 0 || direction > 2) + { + throw std::invalid_argument("SNAP derivative direction must be 0, 1, or 2."); + } + const double dx = direction == 0 ? 1.0 : 0.0; + const double dy = direction == 1 ? 1.0 : 0.0; + const double dz = direction == 2 ? 1.0 : 0.0; + const Dual x(neighbor.x, dx); + const Dual y(neighbor.y, dy); + const Dual z(neighbor.z, dz); + const Dual rsq = x * x + y * y + z * z; + const Dual distance = dual_sqrt(rsq); + const double cutoff_radius = neighbor.cutoff_radius > 0.0 ? neighbor.cutoff_radius : parameters_.rcut; + + UWithDerivative result; + result.derivatives.assign(u_size_, std::complex(0.0, 0.0)); + if (distance.value <= 1.0e-14 || distance.value >= cutoff_radius) + { + return result; + } + const SwitchValue switch_value = switching(distance.value, neighbor, cutoff_radius); + if (switch_value.value == 0.0 || neighbor.weight == 0.0) + { + return result; + } + if (distance.value <= parameters_.rmin) + { + throw std::domain_error("SNAP neighbor distance must be greater than rmin0."); + } + const double scale_value = switch_value.value * neighbor.weight; + const double scale_derivative = switch_value.derivative * distance.derivative * neighbor.weight; + const Dual scale(scale_value, scale_derivative); + + const Dual theta0 = (distance - parameters_.rmin) * parameters_.rfac0 * ModuleBase::PI + / (cutoff_radius - parameters_.rmin); + const Dual z0 = distance / dual_tan(theta0); + const Dual r0inv = 1.0 / dual_sqrt(rsq + z0 * z0); + const Dual ar = r0inv * z0; + const Dual ai = -r0inv * z; + const Dual br = r0inv * y; + const Dual bi = -r0inv * x; + std::vector > raw(u_size_, std::complex(Dual(0.0), Dual(0.0))); + raw[0] = std::complex(Dual(1.0), Dual(0.0)); + const int dimension = parameters_.twojmax + 1; + for (int j = 1; j <= parameters_.twojmax; ++j) + { + for (int mb = 0; 2 * mb <= j; ++mb) + { + raw[u_index(j, 0, mb)] = std::complex(Dual(0.0), Dual(0.0)); + for (int ma = 0; ma < j; ++ma) + { + const std::size_t current = u_index(j, ma, mb); + const std::size_t previous = u_index(j - 1, ma, mb); + const double first_root = rootpq_[static_cast((j - ma) * dimension + (j - mb))]; + const double second_root = rootpq_[static_cast((ma + 1) * dimension + (j - mb))]; + raw[current] += Dual(first_root) * std::complex(ar, -ai) * raw[previous]; + raw[current + 1] = Dual(-second_root) * std::complex(br, -bi) * raw[previous]; + } + } + for (int mb = 0; 2 * mb <= j; ++mb) + { + for (int ma = 0; ma <= j; ++ma) + { + const std::size_t left = u_index(j, ma, mb); + const std::size_t right = u_index(j, j - ma, j - mb); + const double sign = ((ma - mb) % 2 == 0) ? 1.0 : -1.0; + raw[right] = Dual(sign) * std::conj(raw[left]); + } + } + } + for (std::size_t i = 0; i < raw.size(); ++i) + { + const std::complex value = scale * raw[i]; + result.derivatives[i] = std::complex(value.real().derivative, value.imag().derivative); + } + return result; +} + +} // namespace SNAP +} // namespace ModuleMD diff --git a/source/source_md/module_snap/snap_descriptor.cpp b/source/source_md/module_snap/snap_descriptor.cpp new file mode 100644 index 00000000000..9689f3a091a --- /dev/null +++ b/source/source_md/module_snap/snap_descriptor.cpp @@ -0,0 +1,417 @@ +#include "snap_descriptor.h" + +#include "source_base/constants.h" + +#include +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +DescriptorEvaluator::DescriptorEvaluator(const ModelParameters& parameters) + : parameters_(parameters), cg_table_(parameters.twojmax), u_size_(0) +{ + std::string error; + if (!parameters_.valid(error)) + { + throw std::invalid_argument(error); + } + + u_block_.resize(static_cast(parameters_.twojmax + 1), 0); + for (int j = 0; j <= parameters_.twojmax; ++j) + { + u_block_[static_cast(j)] = u_size_; + u_size_ += static_cast((j + 1) * (j + 1)); + } + compute_rootpq(); + + for (int j1 = 0; j1 <= parameters_.twojmax; ++j1) + { + for (int j2 = 0; j2 <= j1; ++j2) + { + for (int j = j1 - j2; j <= std::min(parameters_.twojmax, j1 + j2); j += 2) + { + if (j >= j1) + { + BispectrumIndex index; + index.j1 = j1; + index.j2 = j2; + index.j = j; + bispectrum_indices_.push_back(index); + } + } + } + } +} + +void DescriptorEvaluator::compute_rootpq() +{ + const int dimension = parameters_.twojmax + 1; + rootpq_.assign(static_cast(dimension * dimension), 0.0); + for (int p = 1; p <= parameters_.twojmax; ++p) + { + for (int q = 1; q <= parameters_.twojmax; ++q) + { + rootpq_[static_cast(p * dimension + q)] = std::sqrt(static_cast(p) / q); + } + } +} + +std::size_t DescriptorEvaluator::u_index(const int j, const int ma, const int mb) const +{ + return u_block_[static_cast(j)] + static_cast(mb * (j + 1) + ma); +} + +DescriptorEvaluator::SwitchValue DescriptorEvaluator::switching(const double distance, + const Neighbor& neighbor, + const double cutoff_radius) const +{ + SwitchValue result; + result.value = 1.0; + result.derivative = 0.0; + if (distance >= cutoff_radius) + { + result.value = 0.0; + return result; + } + if (parameters_.switch_flag && distance > parameters_.rmin) + { + const double angle_scale = ModuleBase::PI / (cutoff_radius - parameters_.rmin); + const double angle = (distance - parameters_.rmin) * angle_scale; + result.value = 0.5 * (std::cos(angle) + 1.0); + result.derivative = -0.5 * std::sin(angle) * angle_scale; + } + if (parameters_.switchinner_flag) + { + const double inner_start = neighbor.sinner - neighbor.dinner; + const double inner_end = neighbor.sinner + neighbor.dinner; + if (distance <= inner_start) + { + result.value = 0.0; + result.derivative = 0.0; + } + else if (distance < inner_end) + { + const double angle_scale = 0.5 * ModuleBase::PI / neighbor.dinner; + const double angle = (distance - inner_start) * angle_scale; + const double inner_value = 0.5 * (1.0 - std::cos(angle)); + const double inner_derivative = 0.5 * std::sin(angle) * angle_scale; + result.derivative = result.derivative * inner_value + result.value * inner_derivative; + result.value *= inner_value; + } + } + return result; +} + +std::vector > DescriptorEvaluator::compute_u(const Neighbor& neighbor) const +{ + const double rsq = neighbor.x * neighbor.x + neighbor.y * neighbor.y + neighbor.z * neighbor.z; + const double distance = std::sqrt(rsq); + const double cutoff_radius = neighbor.cutoff_radius > 0.0 ? neighbor.cutoff_radius : parameters_.rcut; + std::vector > result(u_size_, std::complex(0.0, 0.0)); + if (distance <= 1.0e-14 || switching(distance, neighbor, cutoff_radius).value == 0.0) + { + return result; + } + if (distance <= parameters_.rmin) + { + throw std::domain_error("SNAP neighbor distance must be greater than rmin0."); + } + + const double theta0 = (distance - parameters_.rmin) * parameters_.rfac0 * ModuleBase::PI + / (cutoff_radius - parameters_.rmin); + const double z0 = distance / std::tan(theta0); + const double r0inv = 1.0 / std::sqrt(rsq + z0 * z0); + const double ar = r0inv * z0; + const double ai = -r0inv * neighbor.z; + const double br = r0inv * neighbor.y; + const double bi = -r0inv * neighbor.x; + const int dimension = parameters_.twojmax + 1; + + result[0] = std::complex(1.0, 0.0); + for (int j = 1; j <= parameters_.twojmax; ++j) + { + for (int mb = 0; 2 * mb <= j; ++mb) + { + result[u_index(j, 0, mb)] = std::complex(0.0, 0.0); + for (int ma = 0; ma < j; ++ma) + { + const std::size_t current = u_index(j, ma, mb); + const std::size_t previous = u_index(j - 1, ma, mb); + const double first_root = rootpq_[static_cast((j - ma) * dimension + (j - mb))]; + const double second_root = rootpq_[static_cast((ma + 1) * dimension + (j - mb))]; + result[current] += first_root * std::complex(ar, -ai) * result[previous]; + result[current + 1] = -second_root * std::complex(br, -bi) * result[previous]; + } + } + + for (int mb = 0; 2 * mb <= j; ++mb) + { + for (int ma = 0; ma <= j; ++ma) + { + const std::size_t left = u_index(j, ma, mb); + const std::size_t right = u_index(j, j - ma, j - mb); + const double sign = ((ma - mb) % 2 == 0) ? 1.0 : -1.0; + result[right] = sign * std::conj(result[left]); + } + } + } + return result; +} + +std::vector > DescriptorEvaluator::accumulate_u(const std::vector& neighbors, + const std::size_t central_element) const +{ + const std::size_t channels = channel_count(); + const std::size_t self_channel = parameters_.chem_flag ? central_element : 0; + if (self_channel >= channels) + { + throw std::out_of_range("SNAP central element is outside the descriptor channels."); + } + std::vector > result(channels * u_size_, std::complex(0.0, 0.0)); + for (std::size_t channel = 0; channel < channels; ++channel) + { + if (channel != self_channel && !parameters_.wselfall_flag) + { + continue; + } + for (int j = 0; j <= parameters_.twojmax; ++j) + { + for (int ma = 0; ma <= j; ++ma) + { + result[channel * u_size_ + u_index(j, ma, ma)] = parameters_.wself; + } + } + } + + for (std::size_t index = 0; index < neighbors.size(); ++index) + { + const Neighbor& neighbor = neighbors[index]; + const std::size_t channel = parameters_.chem_flag ? neighbor.element_index : 0; + if (channel >= channels) + { + throw std::out_of_range("SNAP neighbor element is outside the descriptor channels."); + } + const double distance = std::sqrt(neighbor.x * neighbor.x + neighbor.y * neighbor.y + neighbor.z * neighbor.z); + const double cutoff_radius = neighbor.cutoff_radius > 0.0 ? neighbor.cutoff_radius : parameters_.rcut; + const double scale = switching(distance, neighbor, cutoff_radius).value * neighbor.weight; + if (scale == 0.0) + { + continue; + } + const std::vector > one_neighbor = compute_u(neighbor); + const std::size_t offset = channel * u_size_; + for (std::size_t i = 0; i < one_neighbor.size(); ++i) + { + result[offset + i] += scale * one_neighbor[i]; + } + } + return result; +} + +std::vector DescriptorEvaluator::compute_bispectrum(const std::vector >& u) const +{ + const std::size_t channels = channel_count(); + if (u.size() != channels * u_size_) + { + throw std::invalid_argument("SNAP u array has an unexpected size."); + } + const std::size_t bispectrum_size = bispectrum_indices_.size(); + std::vector result(channels * channels * channels * bispectrum_size, 0.0); + const auto compute_one = [this, &u](const std::size_t first_channel, + const std::size_t second_channel, + const std::size_t third_channel, + const std::size_t bispectrum_index) { + const BispectrumIndex& index = bispectrum_indices_[bispectrum_index]; + const std::size_t first_offset = first_channel * u_size_; + const std::size_t second_offset = second_channel * u_size_; + const std::size_t third_offset = third_channel * u_size_; + const auto compute_z = [this, &u, &index, first_offset, second_offset](const int ma, const int mb) { + std::complex z(0.0, 0.0); + const int m = 2 * ma - index.j; + const int mp = 2 * mb - index.j; + for (int ma1 = 0; ma1 <= index.j1; ++ma1) + { + const int m1 = 2 * ma1 - index.j1; + for (int ma2 = 0; ma2 <= index.j2; ++ma2) + { + const int m2 = 2 * ma2 - index.j2; + const double c1 = cg_table_.coefficient(index.j1, m1, index.j2, m2, index.j, m); + if (c1 == 0.0) + { + continue; + } + for (int mb1 = 0; mb1 <= index.j1; ++mb1) + { + const int mp1 = 2 * mb1 - index.j1; + for (int mb2 = 0; mb2 <= index.j2; ++mb2) + { + const int mp2 = 2 * mb2 - index.j2; + const double c2 = cg_table_.coefficient(index.j1, mp1, index.j2, mp2, index.j, mp); + if (c2 != 0.0) + { + z += c1 * c2 + * u[first_offset + u_index(index.j1, ma1, mb1)] + * u[second_offset + u_index(index.j2, ma2, mb2)]; + } + } + } + } + } + return z; + }; + std::complex value(0.0, 0.0); + for (int mb = 0; 2 * mb < index.j; ++mb) + { + for (int ma = 0; ma <= index.j; ++ma) + { + value += std::conj(u[third_offset + u_index(index.j, ma, mb)]) * compute_z(ma, mb); + } + } + if (index.j % 2 == 0) + { + const int mb = index.j / 2; + for (int ma = 0; ma < mb; ++ma) + { + value += std::conj(u[third_offset + u_index(index.j, ma, mb)]) * compute_z(ma, mb); + } + value += 0.5 * std::conj(u[third_offset + u_index(index.j, mb, mb)]) * compute_z(mb, mb); + } + return 2.0 * value.real() / (parameters_.bnorm_flag ? index.j + 1.0 : 1.0); + }; + for (std::size_t first_channel = 0; first_channel < channels; ++first_channel) + { + for (std::size_t second_channel = 0; second_channel < channels; ++second_channel) + { + for (std::size_t third_channel = 0; third_channel < channels; ++third_channel) + { + const std::size_t triple = (first_channel * channels + second_channel) * channels + third_channel; + for (std::size_t ib = 0; ib < bispectrum_size; ++ib) + { + result[triple * bispectrum_size + ib] + = compute_one(first_channel, second_channel, third_channel, ib); + } + } + } + } + return result; +} + +std::vector DescriptorEvaluator::compute_bispectrum_derivative( + const std::vector >& u, + const std::vector >& du) const +{ + const std::size_t channels = channel_count(); + if (u.size() != channels * u_size_ || du.size() != channels * u_size_) + { + throw std::invalid_argument("SNAP u derivative array has an unexpected size."); + } + const std::size_t bispectrum_size = bispectrum_indices_.size(); + std::vector result(channels * channels * channels * bispectrum_size, 0.0); + const auto compute_one = [this, &u, &du](const std::size_t first_channel, + const std::size_t second_channel, + const std::size_t third_channel, + const std::size_t bispectrum_index) { + const BispectrumIndex& index = bispectrum_indices_[bispectrum_index]; + const std::size_t first_offset = first_channel * u_size_; + const std::size_t second_offset = second_channel * u_size_; + const std::size_t third_offset = third_channel * u_size_; + const auto compute_z = [this, &u, &du, &index, first_offset, second_offset](const int ma, + const int mb, + std::complex& z, + std::complex& dz) { + z = std::complex(0.0, 0.0); + dz = std::complex(0.0, 0.0); + const int m = 2 * ma - index.j; + const int mp = 2 * mb - index.j; + for (int ma1 = 0; ma1 <= index.j1; ++ma1) + { + const int m1 = 2 * ma1 - index.j1; + for (int ma2 = 0; ma2 <= index.j2; ++ma2) + { + const int m2 = 2 * ma2 - index.j2; + const double c1 = cg_table_.coefficient(index.j1, m1, index.j2, m2, index.j, m); + if (c1 == 0.0) + { + continue; + } + for (int mb1 = 0; mb1 <= index.j1; ++mb1) + { + const int mp1 = 2 * mb1 - index.j1; + for (int mb2 = 0; mb2 <= index.j2; ++mb2) + { + const int mp2 = 2 * mb2 - index.j2; + const double c2 = cg_table_.coefficient(index.j1, mp1, index.j2, mp2, index.j, mp); + if (c2 == 0.0) + { + continue; + } + const double coefficient = c1 * c2; + const std::complex u1 = u[first_offset + u_index(index.j1, ma1, mb1)]; + const std::complex u2 = u[second_offset + u_index(index.j2, ma2, mb2)]; + const std::complex du1 = du[first_offset + u_index(index.j1, ma1, mb1)]; + const std::complex du2 = du[second_offset + u_index(index.j2, ma2, mb2)]; + z += coefficient * u1 * u2; + dz += coefficient * (du1 * u2 + u1 * du2); + } + } + } + } + }; + std::complex derivative(0.0, 0.0); + for (int mb = 0; 2 * mb < index.j; ++mb) + { + for (int ma = 0; ma <= index.j; ++ma) + { + std::complex z(0.0, 0.0); + std::complex dz(0.0, 0.0); + compute_z(ma, mb, z, dz); + const std::complex du3 = du[third_offset + u_index(index.j, ma, mb)]; + derivative += std::conj(du3) * z + + std::conj(u[third_offset + u_index(index.j, ma, mb)]) * dz; + } + } + if (index.j % 2 == 0) + { + const int mb = index.j / 2; + for (int ma = 0; ma < mb; ++ma) + { + std::complex z(0.0, 0.0); + std::complex dz(0.0, 0.0); + compute_z(ma, mb, z, dz); + const std::complex du3 = du[third_offset + u_index(index.j, ma, mb)]; + derivative += std::conj(du3) * z + + std::conj(u[third_offset + u_index(index.j, ma, mb)]) * dz; + } + std::complex z(0.0, 0.0); + std::complex dz(0.0, 0.0); + compute_z(mb, mb, z, dz); + const std::complex du3 = du[third_offset + u_index(index.j, mb, mb)]; + derivative += 0.5 * (std::conj(du3) * z + + std::conj(u[third_offset + u_index(index.j, mb, mb)]) * dz); + } + return 2.0 * derivative.real() / (parameters_.bnorm_flag ? index.j + 1.0 : 1.0); + }; + for (std::size_t first_channel = 0; first_channel < channels; ++first_channel) + { + for (std::size_t second_channel = 0; second_channel < channels; ++second_channel) + { + for (std::size_t third_channel = 0; third_channel < channels; ++third_channel) + { + const std::size_t triple = (first_channel * channels + second_channel) * channels + third_channel; + for (std::size_t ib = 0; ib < bispectrum_size; ++ib) + { + result[triple * bispectrum_size + ib] + = compute_one(first_channel, second_channel, third_channel, ib); + } + } + } + } + return result; +} + +} // namespace SNAP +} // namespace ModuleMD diff --git a/source/source_md/module_snap/snap_descriptor.h b/source/source_md/module_snap/snap_descriptor.h new file mode 100644 index 00000000000..e3759d4eca5 --- /dev/null +++ b/source/source_md/module_snap/snap_descriptor.h @@ -0,0 +1,92 @@ +#ifndef SNAP_DESCRIPTOR_H +#define SNAP_DESCRIPTOR_H + +#include "snap_cg.h" +#include "snap_model.h" + +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +struct Neighbor +{ + double x = 0.0; + double y = 0.0; + double z = 0.0; + double weight = 1.0; + double cutoff_radius = 0.0; + std::size_t element_index = 0; + double sinner = 0.0; + double dinner = 0.0; +}; + +struct UWithDerivative +{ + std::vector > derivatives; +}; + +class DescriptorEvaluator +{ + public: + explicit DescriptorEvaluator(const ModelParameters& parameters); + + std::vector > compute_u(const Neighbor& neighbor) const; + + UWithDerivative compute_u_with_derivative(const Neighbor& neighbor, int direction) const; + + std::vector > accumulate_u(const std::vector& neighbors, + std::size_t central_element) const; + + std::vector compute_bispectrum(const std::vector >& u) const; + + std::vector compute_bispectrum_derivative(const std::vector >& u, + const std::vector >& du) const; + + std::size_t u_size() const + { + return u_size_; + } + + std::size_t bispectrum_size() const + { + return bispectrum_indices_.size(); + } + + std::size_t channel_count() const + { + return parameters_.chem_flag ? parameters_.elements.size() : 1; + } + + private: + struct SwitchValue + { + double value; + double derivative; + }; + struct BispectrumIndex + { + int j1; + int j2; + int j; + }; + + std::size_t u_index(int j, int ma, int mb) const; + SwitchValue switching(double distance, const Neighbor& neighbor, double cutoff_radius) const; + void compute_rootpq(); + + ModelParameters parameters_; + ClebschGordanTable cg_table_; + std::vector rootpq_; + std::vector u_block_; + std::vector bispectrum_indices_; + std::size_t u_size_; +}; + +} // namespace SNAP +} // namespace ModuleMD + +#endif // SNAP_DESCRIPTOR_H diff --git a/source/source_md/module_snap/snap_model.cpp b/source/source_md/module_snap/snap_model.cpp new file mode 100644 index 00000000000..caf8caa866a --- /dev/null +++ b/source/source_md/module_snap/snap_model.cpp @@ -0,0 +1,466 @@ +#include "snap_model.h" +#include +#include +#include +#include + +namespace +{ +bool next_data_line(std::istream& input, std::string& line) +{ + while (std::getline(input, line)) + { + const std::size_t comment = line.find('#'); + if (comment != std::string::npos) + { + line.erase(comment); + } + std::istringstream words(line); + std::string first; + if (words >> first) + { + return true; + } + } + return false; +} + +template +bool read_parameter_value(std::istringstream& words, T& value, const std::string& keyword, std::string& error) +{ + if (!(words >> value) || words.fail()) + { + error = "Invalid SNAP parameter value for keyword: " + keyword; + return false; + } + return true; +} + +bool validate_inner(const ModuleMD::SNAP::ModelParameters& parameters, std::string& error) +{ + if (!parameters.switchinner_flag) + { + return true; + } + const std::size_t count = parameters.elements.size(); + if (parameters.sinner.size() != count || parameters.dinner.size() != count) + { + error = "SNAP switchinnerflag requires sinner and dinner values for every element."; + return false; + } + for (std::size_t i = 0; i < count; ++i) + { + if (!std::isfinite(parameters.sinner[i]) || !std::isfinite(parameters.dinner[i]) + || parameters.dinner[i] <= 0.0) + { + error = "SNAP sinner must be finite and dinner must be finite and positive."; + return false; + } + } + return true; +} + +bool validate_cutoff(const ModuleMD::SNAP::ModelParameters& parameters, std::string& error) +{ + if (!std::isfinite(parameters.rcut) || parameters.rcut <= 0.0) + { + error = "SNAP cutoff must be finite and positive."; + return false; + } + if (!std::isfinite(parameters.rcutfac) || parameters.rcutfac <= 0.0) + { + error = "SNAP rcutfac must be finite and positive."; + return false; + } + return true; +} + +bool validate_pair_cutoffs(const ModuleMD::SNAP::ModelParameters& parameters, std::string& error) +{ + for (std::size_t i = 0; i < parameters.element_radii.size(); ++i) + { + for (std::size_t j = i; j < parameters.element_radii.size(); ++j) + { + const double pair_cutoff = (parameters.element_radii[i] + parameters.element_radii[j]) + * parameters.rcutfac; + if (!std::isfinite(pair_cutoff) || parameters.rmin >= pair_cutoff) + { + error = "SNAP rmin0 must be smaller than the cutoff for element pair " + + parameters.elements[i] + "-" + parameters.elements[j] + "."; + return false; + } + } + } + return true; +} +} + +namespace ModuleMD +{ +namespace SNAP +{ + +bool ModelParameters::valid(std::string& error) const +{ + if (!validate_cutoff(*this, error)) + { + return false; + } + if (!std::isfinite(rmin) || rmin < 0.0 || rmin >= rcut) + { + error = "SNAP rmin must satisfy 0 <= rmin < rcut."; + return false; + } + if (!std::isfinite(rfac0) || rfac0 <= 0.0 || rfac0 > 1.0) + { + error = "SNAP rfac0 must satisfy 0 < rfac0 <= 1."; + return false; + } + if (twojmax < 0) + { + error = "SNAP twojmax must be non-negative."; + return false; + } + if (!std::isfinite(wself)) + { + error = "SNAP self weight must be finite."; + return false; + } + if (elements.empty() || element_weights.size() != elements.size()) + { + error = "SNAP element names and weights must be non-empty and have the same size."; + return false; + } + if (element_radii.size() != elements.size()) + { + error = "SNAP element names and radii must have the same size."; + return false; + } + if (!validate_inner(*this, error)) + { + return false; + } + for (std::size_t i = 0; i < element_weights.size(); ++i) + { + if (!std::isfinite(element_weights[i]) || !std::isfinite(element_radii[i]) || element_radii[i] <= 0.0) + { + error = "SNAP element weights must be finite."; + return false; + } + } + if (!validate_pair_cutoffs(*this, error)) + { + return false; + } + return true; +} + +bool ModelData::load(const std::string& coefficient_file, + const std::string& parameter_file, + const std::vector& requested_elements, + std::string& error) +{ + error.clear(); + parameters = ModelParameters(); + parameters.rfac0 = 0.99363; + coefficients.clear(); + parameters.elements = requested_elements; + if (parameters.elements.empty()) + { + error = "SNAP requires at least one requested element."; + return false; + } + + std::ifstream coeff_stream(coefficient_file.c_str()); + if (!coeff_stream) + { + error = "Cannot open SNAP coefficient file: " + coefficient_file; + return false; + } + + std::string line; + if (!next_data_line(coeff_stream, line)) + { + error = "SNAP coefficient file is empty."; + return false; + } + std::istringstream header(line); + int block_count = 0; + int coefficient_count = 0; + if (!(header >> block_count >> coefficient_count) || block_count <= 0 || coefficient_count <= 0) + { + error = "Invalid SNAP coefficient file header."; + return false; + } + std::string extra; + if (header >> extra) + { + error = "Unexpected extra value in SNAP coefficient file header."; + return false; + } + + std::vector found(parameters.elements.size(), false); + parameters.element_weights.assign(parameters.elements.size(), 0.0); + parameters.element_radii.assign(parameters.elements.size(), 0.0); + coefficients.assign(parameters.elements.size(), std::vector()); + for (int block = 0; block < block_count; ++block) + { + if (!next_data_line(coeff_stream, line)) + { + error = "SNAP coefficient file ended before all element blocks were read."; + return false; + } + std::istringstream element_line(line); + std::string element; + double radius = 0.0; + double weight = 0.0; + if (!(element_line >> element >> radius >> weight)) + { + error = "Invalid SNAP element line."; + return false; + } + if (element_line >> extra) + { + error = "Unexpected extra value in SNAP element line: " + element; + return false; + } + std::size_t element_index = parameters.elements.size(); + for (std::size_t i = 0; i < parameters.elements.size(); ++i) + { + if (parameters.elements[i] == element) + { + element_index = i; + break; + } + } + std::vector values(static_cast(coefficient_count), 0.0); + for (int coefficient = 0; coefficient < coefficient_count; ++coefficient) + { + if (!next_data_line(coeff_stream, line)) + { + error = "SNAP coefficient file ended before all coefficients were read."; + return false; + } + std::istringstream value_line(line); + if (!(value_line >> values[static_cast(coefficient)])) + { + error = "Invalid SNAP coefficient value."; + return false; + } + if (!std::isfinite(values[static_cast(coefficient)]) || value_line >> extra) + { + error = "Invalid or extra SNAP coefficient value for element: " + element; + return false; + } + } + if (element_index < parameters.elements.size()) + { + if (found[element_index]) + { + error = "Duplicate SNAP element block: " + element; + return false; + } + found[element_index] = true; + parameters.element_radii[element_index] = radius; + parameters.element_weights[element_index] = weight; + coefficients[element_index] = values; + } + } + if (next_data_line(coeff_stream, line)) + { + error = "Unexpected data after the SNAP coefficient blocks."; + return false; + } + for (std::size_t i = 0; i < found.size(); ++i) + { + if (!found[i]) + { + error = "SNAP coefficient file does not contain requested element: " + parameters.elements[i]; + return false; + } + } + + std::ifstream parameter_stream(parameter_file.c_str()); + if (!parameter_stream) + { + error = "Cannot open SNAP parameter file: " + parameter_file; + return false; + } + bool has_rcutfac = false; + bool has_twojmax = false; + bool has_sinner = false; + bool has_dinner = false; + while (next_data_line(parameter_stream, line)) + { + std::istringstream words(line); + std::string keyword; + words >> keyword; + if (keyword == "rcutfac") + { + if (!read_parameter_value(words, parameters.rcutfac, keyword, error)) + { + return false; + } + has_rcutfac = true; + } + else if (keyword == "twojmax") + { + if (!read_parameter_value(words, parameters.twojmax, keyword, error)) + { + return false; + } + has_twojmax = true; + } + else if (keyword == "rfac0") + { + if (!read_parameter_value(words, parameters.rfac0, keyword, error)) + { + return false; + } + } + else if (keyword == "rmin0") + { + if (!read_parameter_value(words, parameters.rmin, keyword, error)) + { + return false; + } + } + else if (keyword == "bnormflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.bnorm_flag = value != 0; + } + else if (keyword == "switchflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.switch_flag = value != 0; + } + else if (keyword == "bzeroflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.bzero_flag = value != 0; + } + else if (keyword == "quadraticflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.quadratic_flag = value != 0; + } + else if (keyword == "wselfallflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.wselfall_flag = value != 0; + } + else if (keyword == "chemflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + parameters.chem_flag = value != 0; + } + else if (keyword == "switchinnerflag") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + if (value != 0 && value != 1) + { + error = "SNAP switchinnerflag must be 0 or 1."; + return false; + } + parameters.switchinner_flag = value == 1; + } + else if (keyword == "sinner" || keyword == "dinner") + { + std::vector& values = keyword == "sinner" ? parameters.sinner : parameters.dinner; + values.resize(parameters.elements.size()); + for (std::size_t i = 0; i < values.size(); ++i) + { + if (!read_parameter_value(words, values[i], keyword, error)) + { + return false; + } + } + if (keyword == "sinner") + { + has_sinner = true; + } + else + { + has_dinner = true; + } + } + else if (keyword == "chunksize" || keyword == "parallelthresh") + { + int value = 0; + if (!read_parameter_value(words, value, keyword, error)) + { + return false; + } + if (value <= 0) + { + error = "SNAP " + keyword + " must be positive."; + return false; + } + } + else + { + error = "Unknown SNAP parameter: " + keyword; + return false; + } + std::string extra; + if (words >> extra) + { + error = "Unexpected extra value for SNAP parameter: " + keyword; + return false; + } + } + if (parameters.switchinner_flag && (!has_sinner || !has_dinner)) + { + error = "SNAP switchinnerflag=1 requires sinner and dinner."; + return false; + } + if (!parameters.switchinner_flag && (has_sinner || has_dinner)) + { + error = "SNAP sinner and dinner require switchinnerflag=1."; + return false; + } + if (!has_rcutfac || !has_twojmax) + { + error = "SNAP parameter file must define rcutfac and twojmax."; + return false; + } + double maximum_cutoff = 0.0; + for (std::size_t i = 0; i < parameters.element_radii.size(); ++i) + { + maximum_cutoff = std::max(maximum_cutoff, 2.0 * parameters.element_radii[i] * parameters.rcutfac); + } + parameters.rcut = maximum_cutoff; + return parameters.valid(error); +} + +} // namespace SNAP +} // namespace ModuleMD diff --git a/source/source_md/module_snap/snap_model.h b/source/source_md/module_snap/snap_model.h new file mode 100644 index 00000000000..29d7c2cbe83 --- /dev/null +++ b/source/source_md/module_snap/snap_model.h @@ -0,0 +1,51 @@ +#ifndef SNAP_MODEL_H +#define SNAP_MODEL_H + +#include +#include + +namespace ModuleMD +{ +namespace SNAP +{ + +struct ModelParameters +{ + double rcut = 0.0; + double rcutfac = 0.0; + double rmin = 0.0; + double rfac0 = 0.0; + int twojmax = 0; + bool bnorm_flag = false; + bool chem_flag = false; + bool switch_flag = true; + bool bzero_flag = true; + bool quadratic_flag = false; + bool wselfall_flag = false; + bool switchinner_flag = false; + double wself = 1.0; + std::vector elements; + std::vector element_weights; + std::vector element_radii; + std::vector sinner; + std::vector dinner; + + bool valid(std::string& error) const; +}; + +struct ModelData +{ + ModelParameters parameters; + std::vector > coefficients; + + bool load(const std::string& coefficient_file, + const std::string& parameter_file, + const std::vector& requested_elements, + std::string& error); + +}; + +} // namespace SNAP +} // namespace ModuleMD + +#endif // SNAP_MODEL_H diff --git a/source/source_md/module_snap/unittest/CMakeLists.txt b/source/source_md/module_snap/unittest/CMakeLists.txt new file mode 100644 index 00000000000..86495fdb4ed --- /dev/null +++ b/source/source_md/module_snap/unittest/CMakeLists.txt @@ -0,0 +1,21 @@ +AddTest( + TARGET MODULE_MD_SNAP_model + SOURCES test_snap_model.cpp ../snap_model.cpp +) + +AddTest( + TARGET MODULE_MD_SNAP_cg + SOURCES test_snap_cg.cpp ../snap_cg.cpp +) + +AddTest( + TARGET MODULE_MD_SNAP_descriptor + LIBS base device + SOURCES test_snap_descriptor.cpp ../snap_descriptor.cpp ../snap_derivative.cpp ../snap_model.cpp ../snap_cg.cpp +) + +AddTest( + TARGET MODULE_MD_SNAP_calculator + LIBS base device + SOURCES test_snap_calculator.cpp ../snap_calculator.cpp ../snap_descriptor.cpp ../snap_derivative.cpp ../snap_model.cpp ../snap_cg.cpp +) diff --git a/source/source_md/module_snap/unittest/test_snap_calculator.cpp b/source/source_md/module_snap/unittest/test_snap_calculator.cpp new file mode 100644 index 00000000000..70dd85dae55 --- /dev/null +++ b/source/source_md/module_snap/unittest/test_snap_calculator.cpp @@ -0,0 +1,460 @@ +#include "../snap_calculator.h" + +#include + +#include +#include + +namespace +{ + +ModuleMD::SNAP::ModelData make_model() +{ + ModuleMD::SNAP::ModelData model; + model.parameters.rcut = 4.5; + model.parameters.rcutfac = 1.0; + model.parameters.rmin = 0.0; + model.parameters.rfac0 = 0.99363; + model.parameters.twojmax = 0; + model.parameters.elements.push_back("Si"); + model.parameters.element_weights.push_back(1.0); + model.parameters.element_radii.push_back(2.25); + model.coefficients.push_back(std::vector{2.0, 3.0}); + return model; +} + +ModuleMD::SNAP::ModelData make_multi_element_model() +{ + ModuleMD::SNAP::ModelData model; + model.parameters.rcut = 4.5; + model.parameters.rcutfac = 1.0; + model.parameters.rmin = 0.0; + model.parameters.rfac0 = 0.99363; + model.parameters.twojmax = 0; + model.parameters.bzero_flag = false; + model.parameters.elements.push_back("Si"); + model.parameters.elements.push_back("C"); + model.parameters.element_weights.push_back(1.0); + model.parameters.element_weights.push_back(0.8); + model.parameters.element_radii.push_back(2.25); + model.parameters.element_radii.push_back(1.7); + model.coefficients.push_back(std::vector{2.0, 3.0}); + model.coefficients.push_back(std::vector{5.0, 7.0}); + return model; +} + +ModuleMD::SNAP::ModelData make_chemical_model() +{ + ModuleMD::SNAP::ModelData model; + model.parameters.rcut = 4.5; + model.parameters.rcutfac = 1.0; + model.parameters.rmin = 0.0; + model.parameters.rfac0 = 0.99363; + model.parameters.twojmax = 0; + model.parameters.bnorm_flag = true; + model.parameters.chem_flag = true; + model.parameters.wselfall_flag = true; + model.parameters.elements.push_back("Si"); + model.parameters.elements.push_back("C"); + model.parameters.element_weights.push_back(1.0); + model.parameters.element_weights.push_back(0.8); + model.parameters.element_radii.push_back(2.25); + model.parameters.element_radii.push_back(1.7); + model.coefficients.push_back(std::vector(9, 0.0)); + model.coefficients.push_back(std::vector(9, 0.0)); + model.coefficients[0][1] = 2.0; + model.coefficients[1][1] = 3.0; + return model; +} + +ModuleMD::SNAP::ModelData make_angular_quadratic_model() +{ + ModuleMD::SNAP::ModelData model = make_model(); + model.parameters.twojmax = 2; + model.parameters.quadratic_flag = true; + model.coefficients[0] = {0.2, 0.3, -0.1, 0.5, 0.2, -0.4, + 0.25, 0.0, 0.0, 0.0, -0.15, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0}; + return model; +} + +struct SmallCellResult +{ + double energy = 0.0; + double virial_xx = 0.0; + std::array, 3> forces = {}; +}; + +SmallCellResult evaluate_small_cell(const ModuleMD::SNAP::Calculator& calculator, + const std::array, 3>& positions) +{ + SmallCellResult total; + ModuleMD::SNAP::CellAccumulator accumulator(positions.size()); + for (std::size_t center = 0; center < positions.size(); ++center) + { + std::vector neighbors; + std::vector neighbor_indices; + for (std::size_t atom = 0; atom < positions.size(); ++atom) + { + if (atom == center) + { + continue; + } + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = positions[atom][0] - positions[center][0]; + neighbor.y = positions[atom][1] - positions[center][1]; + neighbor.z = positions[atom][2] - positions[center][2]; + neighbor.cutoff_radius = 4.5; + neighbors.push_back(neighbor); + neighbor_indices.push_back(atom); + } + const ModuleMD::SNAP::ForceResult result = calculator.compute_with_analytic_force(neighbors, 0); + accumulator.add(center, neighbor_indices, neighbors, result); + } + total.energy = accumulator.energy(); + total.virial_xx = accumulator.virial()[0]; + for (std::size_t atom = 0; atom < positions.size(); ++atom) + { + total.forces[atom] = accumulator.forces()[atom]; + } + return total; +} + +TEST(SnapCalculator, ComputesLinearAtomicEnergy) +{ + const ModuleMD::SNAP::Calculator calculator(make_model()); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 1.0; + const ModuleMD::SNAP::EnergyResult result = + calculator.compute(std::vector(1, neighbor), 0); + + ASSERT_EQ(result.descriptor.size(), static_cast(1)); + EXPECT_TRUE(std::isfinite(result.descriptor[0])); + EXPECT_DOUBLE_EQ(result.energy, 2.0 + 3.0 * result.descriptor[0]); +} + +TEST(SnapCalculator, SupportsMultipleCentralElements) +{ + const ModuleMD::SNAP::Calculator calculator(make_multi_element_model()); + const std::vector neighbors; + const ModuleMD::SNAP::EnergyResult silicon = calculator.compute(neighbors, 0); + const ModuleMD::SNAP::EnergyResult carbon = calculator.compute(neighbors, 1); + + EXPECT_DOUBLE_EQ(silicon.energy, 2.0 + 3.0 * silicon.descriptor[0]); + EXPECT_DOUBLE_EQ(carbon.energy, 5.0 + 7.0 * carbon.descriptor[0]); + EXPECT_NE(silicon.energy, carbon.energy); +} + +TEST(SnapCalculator, SupportsChemicalChannels) +{ + const ModuleMD::SNAP::Calculator calculator(make_chemical_model()); + ModuleMD::SNAP::Neighbor silicon; + silicon.x = 1.0; + silicon.cutoff_radius = 4.5; + silicon.element_index = 0; + ModuleMD::SNAP::Neighbor carbon = silicon; + carbon.y = 1.0; + carbon.element_index = 1; + const std::vector neighbors = {silicon, carbon}; + + const ModuleMD::SNAP::ForceResult analytic = calculator.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult finite_difference + = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + + EXPECT_EQ(analytic.energy.descriptor.size(), static_cast(8)); + EXPECT_NEAR(analytic.energy.energy, finite_difference.energy.energy, 1.0e-12); + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_NEAR(analytic.center_force[static_cast(direction)], + finite_difference.center_force[static_cast(direction)], + 1.0e-5); + EXPECT_NEAR(analytic.neighbor_forces[0][static_cast(direction)], + finite_difference.neighbor_forces[0][static_cast(direction)], + 1.0e-5); + } +} + +TEST(SnapCalculator, ChemicalAngularChannelsFollowElementOrder) +{ + ModuleMD::SNAP::ModelData model = make_chemical_model(); + model.parameters.twojmax = 2; + const std::vector zero_coefficients(41, 0.0); + model.coefficients.assign(2, zero_coefficients); + // Bispectrum index 2 couples the Si and C angular channels. + const std::size_t original_index = 1 + 2 * 5 + 2; + model.coefficients[0][original_index] = 1.7; + const ModuleMD::SNAP::Calculator calculator(model); + + ModuleMD::SNAP::Neighbor silicon; + silicon.x = 1.0; + silicon.y = 0.1; + silicon.z = 0.2; + silicon.cutoff_radius = 4.5; + silicon.element_index = 0; + ModuleMD::SNAP::Neighbor carbon = silicon; + carbon.x = 0.3; + carbon.y = 1.1; + carbon.z = -0.4; + carbon.weight = 0.8; + carbon.element_index = 1; + const std::vector neighbors = {silicon, carbon}; + const ModuleMD::SNAP::ForceResult original = calculator.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult finite = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + ModuleMD::SNAP::Neighbor rotated_carbon = carbon; + rotated_carbon.x = carbon.y; + rotated_carbon.y = carbon.x; + const std::vector rotated_neighbors = {silicon, rotated_carbon}; + const double rotated_energy = calculator.compute(rotated_neighbors, 0).energy; + EXPECT_GT(std::abs(original.energy.energy - rotated_energy), 1.0e-6); + + ModuleMD::SNAP::ModelData reordered = model; + reordered.parameters.elements = {"C", "Si"}; + reordered.parameters.element_weights = {0.8, 1.0}; + reordered.parameters.element_radii = {1.7, 2.25}; + reordered.coefficients.assign(2, zero_coefficients); + // The same physical channel triple is (Si, C, Si) = (1, 0, 1) after reordering. + const std::size_t reordered_index = 1 + 5 * 5 + 2; + reordered.coefficients[1][reordered_index] = 1.7; + const ModuleMD::SNAP::Calculator reordered_calculator(reordered); + silicon.element_index = 1; + carbon.element_index = 0; + const std::vector reordered_neighbors = {silicon, carbon}; + const ModuleMD::SNAP::ForceResult permuted + = reordered_calculator.compute_with_analytic_force(reordered_neighbors, 1); + + EXPECT_NEAR(original.energy.energy, permuted.energy.energy, 1.0e-10); + for (std::size_t atom = 0; atom < neighbors.size(); ++atom) + { + double force_size = 0.0; + for (int direction = 0; direction < 3; ++direction) + { + const std::size_t component = static_cast(direction); + const double force = original.neighbor_forces[atom][component]; + force_size += std::abs(force); + EXPECT_NEAR(force, permuted.neighbor_forces[atom][component], 1.0e-10); + EXPECT_NEAR(force, finite.neighbor_forces[atom][component], 1.0e-5); + } + EXPECT_GT(force_size, 1.0e-6); + } +} + +TEST(SnapCalculator, RejectsMismatchedCoefficientCount) +{ + ModuleMD::SNAP::ModelData model = make_model(); + model.coefficients[0].push_back(4.0); + try + { + const ModuleMD::SNAP::Calculator calculator(model); + FAIL() << "Expected a coefficient-count error"; + } + catch (const std::invalid_argument& error) + { + EXPECT_NE(std::string(error.what()).find("expected 2, found 3"), std::string::npos); + } +} + +TEST(SnapCalculator, CachedSelfDescriptorMatchesEnergyAndForce) +{ + ModuleMD::SNAP::ModelData model = make_model(); + const ModuleMD::SNAP::Calculator calculator(model); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 1.0; + neighbor.y = 0.4; + const std::vector neighbors(1, neighbor); + const ModuleMD::SNAP::EnergyResult energy = calculator.compute(neighbors, 0); + const ModuleMD::SNAP::ForceResult analytic = calculator.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult finite = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + EXPECT_DOUBLE_EQ(energy.energy, analytic.energy.energy); + ASSERT_EQ(energy.descriptor.size(), analytic.energy.descriptor.size()); + EXPECT_DOUBLE_EQ(energy.descriptor[0], analytic.energy.descriptor[0]); + for (int direction = 0; direction < 3; ++direction) + { + const std::size_t index = static_cast(direction); + EXPECT_NEAR(analytic.neighbor_forces[0][index], finite.neighbor_forces[0][index], 1.0e-6); + } +} + +TEST(SnapCalculator, RejectsInvalidCentralElement) +{ + const ModuleMD::SNAP::Calculator calculator(make_model()); + EXPECT_THROW(calculator.compute(std::vector(), 1), std::out_of_range); +} + +TEST(SnapCalculator, FiniteDifferenceForcesConservePairContribution) +{ + const ModuleMD::SNAP::Calculator calculator(make_model()); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 1.0; + const ModuleMD::SNAP::ForceResult result = calculator.compute_with_finite_difference_force( + std::vector(1, neighbor), 0, 1.0e-5); + + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_DOUBLE_EQ(result.center_force[static_cast(direction)] + + result.neighbor_forces[0][static_cast(direction)], + 0.0); + } +} + +TEST(SnapCalculator, AnalyticForceMatchesFiniteDifference) +{ + const ModuleMD::SNAP::Calculator calculator(make_model()); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 1.0; + neighbor.y = 0.4; + neighbor.z = -0.2; + const std::vector neighbors(1, neighbor); + const ModuleMD::SNAP::ForceResult analytic = calculator.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult finite = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_NEAR(analytic.neighbor_forces[0][static_cast(direction)], + finite.neighbor_forces[0][static_cast(direction)], + 1.0e-6); + } +} + +TEST(SnapCalculator, InnerSwitchMatchesFiniteDifference) +{ + ModuleMD::SNAP::ModelData model = make_model(); + model.parameters.twojmax = 2; + model.parameters.switchinner_flag = true; + model.parameters.sinner.push_back(1.5); + model.parameters.dinner.push_back(0.4); + model.coefficients[0] = {0.2, 0.3, -0.1, 0.5, 0.2, -0.4}; + const ModuleMD::SNAP::Calculator calculator(model); + ModuleMD::SNAP::Neighbor first; + first.x = 1.4; + first.y = 0.3; + first.z = 0.2; + first.sinner = 1.5; + first.dinner = 0.4; + ModuleMD::SNAP::Neighbor second = first; + second.x = -0.3; + second.y = 1.7; + second.z = -0.2; + const std::vector neighbors = {first, second}; + const ModuleMD::SNAP::ForceResult analytic = calculator.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult finite = calculator.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + for (std::size_t i = 0; i < neighbors.size(); ++i) + { + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_NEAR(analytic.neighbor_forces[i][static_cast(direction)], + finite.neighbor_forces[i][static_cast(direction)], 1.0e-5); + } + } + + model.parameters.quadratic_flag = true; + model.coefficients[0].resize(21, 0.0); + model.coefficients[0][6] = 0.25; + model.coefficients[0][10] = -0.15; + const ModuleMD::SNAP::Calculator quadratic(model); + const ModuleMD::SNAP::ForceResult quadratic_analytic = quadratic.compute_with_analytic_force(neighbors, 0); + const ModuleMD::SNAP::ForceResult quadratic_finite = quadratic.compute_with_finite_difference_force(neighbors, 0, 1.0e-6); + for (std::size_t i = 0; i < neighbors.size(); ++i) + { + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_NEAR(quadratic_analytic.neighbor_forces[i][static_cast(direction)], + quadratic_finite.neighbor_forces[i][static_cast(direction)], 1.0e-5); + } + } +} + +TEST(SnapCalculator, InnerSwitchSkipsSingularRadius) +{ + ModuleMD::SNAP::ModelData model = make_model(); + model.parameters.rmin = 0.5; + model.parameters.switchinner_flag = true; + model.parameters.sinner.push_back(1.5); + model.parameters.dinner.push_back(0.4); + const ModuleMD::SNAP::Calculator calculator(model); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 0.5; + neighbor.sinner = 1.5; + neighbor.dinner = 0.4; + const std::vector neighbors(1, neighbor); + EXPECT_NO_THROW(calculator.compute_with_analytic_force(neighbors, 0)); + model.parameters.switchinner_flag = false; + const ModuleMD::SNAP::Calculator bare(model); + EXPECT_THROW(bare.compute(neighbors, 0), std::domain_error); + EXPECT_THROW(bare.compute_with_analytic_force(neighbors, 0), std::domain_error); + + model.parameters.twojmax = 2; + model.coefficients[0] = {0.2, 0.3, -0.1, 0.5, 0.2, -0.4}; + const ModuleMD::SNAP::Calculator angular(model); + neighbor.x = 0.55; + neighbor.y = 0.12; + const std::vector near_neighbors(1, neighbor); + const ModuleMD::SNAP::ForceResult near_analytic = angular.compute_with_analytic_force(near_neighbors, 0); + const ModuleMD::SNAP::ForceResult near_finite = angular.compute_with_finite_difference_force(near_neighbors, 0, 1.0e-6); + EXPECT_TRUE(std::isfinite(near_analytic.energy.energy)); + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_TRUE(std::isfinite(near_analytic.neighbor_forces[0][static_cast(direction)])); + EXPECT_NEAR(near_analytic.neighbor_forces[0][static_cast(direction)], + near_finite.neighbor_forces[0][static_cast(direction)], 1.0e-4); + } +} + +TEST(SnapCalculator, ZeroWeightSkipsInvalidRadius) +{ + ModuleMD::SNAP::ModelData model = make_model(); + model.parameters.rmin = 0.5; + const ModuleMD::SNAP::Calculator calculator(model); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 0.4; + neighbor.weight = 0.0; + const std::vector neighbors(1, neighbor); + const std::vector empty; + const double isolated_energy = calculator.compute(empty, 0).energy; + const ModuleMD::SNAP::ForceResult result = calculator.compute_with_analytic_force(neighbors, 0); + EXPECT_DOUBLE_EQ(result.energy.energy, isolated_energy); + for (int direction = 0; direction < 3; ++direction) + { + EXPECT_DOUBLE_EQ(result.center_force[static_cast(direction)], 0.0); + EXPECT_DOUBLE_EQ(result.neighbor_forces[0][static_cast(direction)], 0.0); + } + neighbor.weight = 1.0; + const std::vector active(1, neighbor); + EXPECT_THROW(calculator.compute(active, 0), std::domain_error); + EXPECT_THROW(calculator.compute_with_analytic_force(active, 0), std::domain_error); +} + +TEST(SnapCalculator, AngularQuadraticForcesMatchLammpsAndStrain) +{ + const ModuleMD::SNAP::Calculator calculator(make_angular_quadratic_model()); + const std::array, 3> positions = {{{0.0, 0.0, 0.0}, + {1.058354, 0.0, 0.0}, + {0.2645885, 1.3229425, 0.0}}}; + const SmallCellResult result = evaluate_small_cell(calculator, positions); + const double lammps_forces[3][2] = {{67.2928038408958, 57.5476436778783}, + {-92.1399449575303, 60.5944497536836}, + {24.8471411166345, -118.142093431562}}; + EXPECT_NEAR(result.energy, 54.909561, 1.0e-5); + for (std::size_t atom = 0; atom < positions.size(); ++atom) + { + for (int direction = 0; direction < 2; ++direction) + { + EXPECT_NEAR(result.forces[atom][static_cast(direction)], + lammps_forces[atom][direction], 2.0e-4); + } + } + + const double strain = 1.0e-6; + std::array, 3> expanded = positions; + std::array, 3> contracted = positions; + for (std::size_t atom = 0; atom < positions.size(); ++atom) + { + expanded[atom][0] *= 1.0 + strain; + contracted[atom][0] *= 1.0 - strain; + } + const double expanded_energy = evaluate_small_cell(calculator, expanded).energy; + const double contracted_energy = evaluate_small_cell(calculator, contracted).energy; + const double strain_derivative = (expanded_energy - contracted_energy) / (2.0 * strain); + EXPECT_NEAR(result.virial_xx, -strain_derivative, 1.0e-3); +} + +} // namespace diff --git a/source/source_md/module_snap/unittest/test_snap_cg.cpp b/source/source_md/module_snap/unittest/test_snap_cg.cpp new file mode 100644 index 00000000000..022dda681d7 --- /dev/null +++ b/source/source_md/module_snap/unittest/test_snap_cg.cpp @@ -0,0 +1,31 @@ +#include "../snap_cg.h" + +#include + +#include + +namespace +{ + +TEST(SnapCG, HighestWeightCoefficientIsOne) +{ + const ModuleMD::SNAP::ClebschGordanTable table(2); + EXPECT_NEAR(table.coefficient(1, 1, 1, 1, 2, 2), 1.0, 1.0e-14); +} + +TEST(SnapCG, SingletHasExpectedSigns) +{ + const ModuleMD::SNAP::ClebschGordanTable table(2); + const double value = 1.0 / std::sqrt(2.0); + EXPECT_NEAR(table.coefficient(1, 1, 1, -1, 0, 0), value, 1.0e-14); + EXPECT_NEAR(table.coefficient(1, -1, 1, 1, 0, 0), -value, 1.0e-14); +} + +TEST(SnapCG, SelectionRulesReturnZero) +{ + const ModuleMD::SNAP::ClebschGordanTable table(2); + EXPECT_DOUBLE_EQ(table.coefficient(1, 1, 1, 1, 0, 0), 0.0); + EXPECT_DOUBLE_EQ(table.coefficient(1, 1, 1, -1, 2, 1), 0.0); +} + +} // namespace diff --git a/source/source_md/module_snap/unittest/test_snap_descriptor.cpp b/source/source_md/module_snap/unittest/test_snap_descriptor.cpp new file mode 100644 index 00000000000..176132000a7 --- /dev/null +++ b/source/source_md/module_snap/unittest/test_snap_descriptor.cpp @@ -0,0 +1,106 @@ +#include "../snap_descriptor.h" + +#include + +#include + +namespace +{ + +ModuleMD::SNAP::ModelParameters make_parameters() +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 4.5; + parameters.rcutfac = 1.0; + parameters.rmin = 0.0; + parameters.rfac0 = 0.99363; + parameters.twojmax = 2; + parameters.bzero_flag = false; + parameters.elements.push_back("Si"); + parameters.element_weights.push_back(1.0); + parameters.element_radii.push_back(2.25); + return parameters; +} + +TEST(SnapDescriptor, AllocatesExpectedArraySizes) +{ + const ModuleMD::SNAP::DescriptorEvaluator evaluator(make_parameters()); + EXPECT_EQ(evaluator.u_size(), static_cast(1 * 1 + 2 * 2 + 3 * 3)); + EXPECT_EQ(evaluator.bispectrum_size(), static_cast(5)); +} + +TEST(SnapDescriptor, SelfContributionIsDiagonal) +{ + const ModuleMD::SNAP::DescriptorEvaluator evaluator(make_parameters()); + const std::vector neighbors; + const std::vector > u = evaluator.accumulate_u(neighbors, 0); + + EXPECT_DOUBLE_EQ(u[0].real(), 1.0); + EXPECT_DOUBLE_EQ(u[0].imag(), 0.0); + EXPECT_DOUBLE_EQ(u[4].real(), 1.0); + EXPECT_DOUBLE_EQ(u[4].imag(), 0.0); + EXPECT_DOUBLE_EQ(u[13].real(), 1.0); + EXPECT_DOUBLE_EQ(u[13].imag(), 0.0); +} + +TEST(SnapDescriptor, OutsideCutoffDoesNotChangeEnvironment) +{ + const ModuleMD::SNAP::DescriptorEvaluator evaluator(make_parameters()); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 5.0; + const std::vector neighbors(1, neighbor); + const std::vector > u = evaluator.accumulate_u(neighbors, 0); + const std::vector no_neighbors; + const std::vector > self = evaluator.accumulate_u(no_neighbors, 0); + + ASSERT_EQ(u.size(), self.size()); + for (std::size_t i = 0; i < u.size(); ++i) + { + EXPECT_NEAR(u[i].real(), self[i].real(), 1.0e-14); + EXPECT_NEAR(u[i].imag(), self[i].imag(), 1.0e-14); + } +} + +TEST(SnapDescriptor, BispectrumIsFinite) +{ + const ModuleMD::SNAP::DescriptorEvaluator evaluator(make_parameters()); + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = 1.2; + neighbor.y = 0.4; + neighbor.z = -0.8; + const std::vector neighbors(1, neighbor); + const std::vector > u = evaluator.accumulate_u(neighbors, 0); + const std::vector bispectrum = evaluator.compute_bispectrum(u); + + for (std::size_t i = 0; i < bispectrum.size(); ++i) + { + EXPECT_TRUE(std::isfinite(bispectrum[i])); + } +} + +TEST(SnapDescriptor, MatchesLammpsSingleElementReference) +{ + const ModuleMD::SNAP::DescriptorEvaluator evaluator(make_parameters()); + const double positions[3][3] = {{1.0, 1.0, 1.0}, {2.0, 1.0, 1.0}, {1.0, 3.0, 1.0}}; + std::vector neighbors; + for (int atom = 1; atom < 3; ++atom) + { + ModuleMD::SNAP::Neighbor neighbor; + neighbor.x = positions[atom][0] - positions[0][0]; + neighbor.y = positions[atom][1] - positions[0][1]; + neighbor.z = positions[atom][2] - positions[0][2]; + neighbor.cutoff_radius = 4.0; + neighbors.push_back(neighbor); + } + const std::vector > u = evaluator.accumulate_u(neighbors, 0); + const std::vector bispectrum = evaluator.compute_bispectrum(u); + const double reference[5] = { + 13.0368348609, 15.0992648922, 11.7631725456, 13.7069739819, 6.07477849248}; + ASSERT_EQ(bispectrum.size(), static_cast(5)); + for (int index = 0; index < 5; ++index) + { + EXPECT_NEAR(bispectrum[static_cast(index)], reference[index], 1.0e-9); + } +} + +} // namespace diff --git a/source/source_md/module_snap/unittest/test_snap_model.cpp b/source/source_md/module_snap/unittest/test_snap_model.cpp new file mode 100644 index 00000000000..118f4fa6a9d --- /dev/null +++ b/source/source_md/module_snap/unittest/test_snap_model.cpp @@ -0,0 +1,295 @@ +#include "../snap_model.h" + +#include + +#include +#include + +namespace +{ + +TEST(SnapModel, AcceptsCompleteParameters) +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 4.5; + parameters.rcutfac = 1.0; + parameters.rmin = 0.0; + parameters.rfac0 = 0.99363; + parameters.twojmax = 8; + parameters.elements.push_back("Si"); + parameters.element_weights.push_back(1.0); + parameters.element_radii.push_back(2.25); + + std::string error; + EXPECT_TRUE(parameters.valid(error)); + EXPECT_TRUE(error.empty()); +} + +TEST(SnapModel, RejectsInconsistentElementData) +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 4.5; + parameters.rcutfac = 1.0; + parameters.rmin = 0.0; + parameters.rfac0 = 0.99363; + parameters.elements.push_back("Si"); + parameters.element_radii.push_back(2.25); + + std::string error; + EXPECT_FALSE(parameters.valid(error)); + EXPECT_FALSE(error.empty()); +} + +TEST(SnapModel, RejectsInvalidRadialRange) +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 3.0; + parameters.rcutfac = 1.0; + parameters.rmin = 3.0; + parameters.rfac0 = 1.0; + parameters.elements.push_back("Si"); + parameters.element_weights.push_back(1.0); + parameters.element_radii.push_back(2.25); + + std::string error; + EXPECT_FALSE(parameters.valid(error)); + EXPECT_FALSE(error.empty()); +} + +TEST(SnapModel, ValidatesRfac0Range) +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 4.5; + parameters.rcutfac = 1.0; + parameters.elements.push_back("Si"); + parameters.element_weights.push_back(1.0); + parameters.element_radii.push_back(2.25); + std::string error; + + parameters.rfac0 = 1.0; + EXPECT_TRUE(parameters.valid(error)) << error; + parameters.rfac0 = 0.0; + EXPECT_FALSE(parameters.valid(error)); + EXPECT_NE(error.find("0 < rfac0 <= 1"), std::string::npos); + parameters.rfac0 = 1.2; + EXPECT_FALSE(parameters.valid(error)); + EXPECT_NE(error.find("0 < rfac0 <= 1"), std::string::npos); +} + +TEST(SnapModel, RejectsRminOutsideSmallPairCutoff) +{ + ModuleMD::SNAP::ModelParameters parameters; + parameters.rcut = 4.0; + parameters.rcutfac = 1.0; + parameters.rmin = 1.5; + parameters.rfac0 = 0.99363; + parameters.elements = {"Small", "Large"}; + parameters.element_weights = {1.0, 1.0}; + parameters.element_radii = {0.5, 2.0}; + std::string error; + EXPECT_FALSE(parameters.valid(error)); + EXPECT_NE(error.find("Small-Small"), std::string::npos); +} + +TEST(SnapModel, LoadsLammpsStyleCoefficientAndParameterFiles) +{ + const std::string coefficient_file = testing::TempDir() + "snap_coeff_test.txt"; + const std::string parameter_file = testing::TempDir() + "snap_param_test.txt"; + { + std::ofstream output(coefficient_file.c_str()); + output << "1 3\n"; + output << "Si 2.25 1.0\n"; + output << "-0.5\n1.25\n2.5\n"; + } + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.2\n"; + output << "twojmax 4\n"; + output << "rfac0 0.99363\n"; + output << "rmin0 0.0\n"; + output << "bnormflag 0\n"; + output << "chemflag 0\n"; + } + + ModuleMD::SNAP::ModelData model; + std::string error; + EXPECT_TRUE(model.load(coefficient_file, parameter_file, std::vector(1, "Si"), error)); + EXPECT_TRUE(error.empty()); + ASSERT_EQ(model.coefficients[0].size(), static_cast(3)); + EXPECT_DOUBLE_EQ(model.parameters.element_weights[0], 1.0); + EXPECT_DOUBLE_EQ(model.parameters.element_radii[0], 2.25); + EXPECT_DOUBLE_EQ(model.parameters.rcut, 5.4); + EXPECT_DOUBLE_EQ(model.coefficients[0][1], 1.25); + + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.2\ntwojmax 4\nrfac0 1.2\n"; + } + EXPECT_FALSE(model.load(coefficient_file, parameter_file, std::vector(1, "Si"), error)); + EXPECT_NE(error.find("0 < rfac0 <= 1"), std::string::npos); + + std::remove(coefficient_file.c_str()); + std::remove(parameter_file.c_str()); +} + +TEST(SnapModel, LoadsMultipleElementBlocks) +{ + const std::string coefficient_file = testing::TempDir() + "snap_coeff_multi_test.txt"; + const std::string parameter_file = testing::TempDir() + "snap_param_multi_test.txt"; + { + std::ofstream output(coefficient_file.c_str()); + output << "2 2\n"; + output << "Si 2.25 1.0\n"; + output << "-0.5\n1.25\n"; + output << "C 1.7 0.8\n"; + output << "-0.25\n2.5\n"; + } + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\n"; + output << "twojmax 0\n"; + output << "rfac0 0.99363\n"; + output << "rmin0 0.0\n"; + output << "chemflag 0\n"; + } + + ModuleMD::SNAP::ModelData model; + std::string error; + const std::vector elements = {"Si", "C"}; + EXPECT_TRUE(model.load(coefficient_file, parameter_file, elements, error)); + EXPECT_TRUE(error.empty()); + ASSERT_EQ(model.coefficients.size(), static_cast(2)); + EXPECT_DOUBLE_EQ(model.parameters.element_weights[1], 0.8); + EXPECT_DOUBLE_EQ(model.parameters.element_radii[1], 1.7); + EXPECT_DOUBLE_EQ(model.coefficients[1][1], 2.5); + + std::remove(coefficient_file.c_str()); + std::remove(parameter_file.c_str()); +} + +TEST(SnapModel, LoadsChemicalCoefficientBlocks) +{ + const std::string coefficient_file = testing::TempDir() + "snap_coeff_chemical_test.txt"; + const std::string parameter_file = testing::TempDir() + "snap_param_chemical_test.txt"; + { + std::ofstream output(coefficient_file.c_str()); + output << "2 9\n"; + output << "Si 2.25 1.0\n"; + for (int i = 0; i < 9; ++i) + { + output << (i == 1 ? 1.0 : 0.0) << "\n"; + } + output << "C 1.7 0.8\n"; + for (int i = 0; i < 9; ++i) + { + output << (i == 1 ? 2.0 : 0.0) << "\n"; + } + } + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\n"; + output << "twojmax 0\n"; + output << "rfac0 0.99363\n"; + output << "rmin0 0.0\n"; + output << "chemflag 1\n"; + output << "bnormflag 1\n"; + output << "wselfallflag 1\n"; + } + + ModuleMD::SNAP::ModelData model; + std::string error; + const std::vector elements = {"Si", "C"}; + EXPECT_TRUE(model.load(coefficient_file, parameter_file, elements, error)); + EXPECT_TRUE(error.empty()); + ASSERT_EQ(model.coefficients.size(), static_cast(2)); + EXPECT_TRUE(model.parameters.chem_flag); + EXPECT_EQ(model.coefficients[0].size(), static_cast(9)); + EXPECT_DOUBLE_EQ(model.coefficients[1][1], 2.0); + + const std::vector reversed_elements = {"C", "Si"}; + ModuleMD::SNAP::ModelData reversed; + ASSERT_TRUE(reversed.load(coefficient_file, parameter_file, reversed_elements, error)) << error; + EXPECT_EQ(reversed.parameters.elements[0], "C"); + EXPECT_EQ(reversed.parameters.elements[1], "Si"); + EXPECT_DOUBLE_EQ(reversed.coefficients[0][1], 2.0); + EXPECT_DOUBLE_EQ(reversed.coefficients[1][1], 1.0); + + std::remove(coefficient_file.c_str()); + std::remove(parameter_file.c_str()); +} + +TEST(SnapModel, ReadsInnerSwitchAndRejectsUnknownParameter) +{ + const std::string coefficient_file = testing::TempDir() + "snap_inner_coeff.txt"; + const std::string parameter_file = testing::TempDir() + "snap_inner_param.txt"; + { + std::ofstream output(coefficient_file.c_str()); + output << "1 2\nSi 2.0 1.0\n0.0\n1.0\n"; + } + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\ntwojmax 0\nswitchinnerflag 1\n"; + output << "sinner 1.5\ndinner 0.4\nchunksize 1024\n"; + } + ModuleMD::SNAP::ModelData model; + std::string error; + const std::vector elements(1, "Si"); + ASSERT_TRUE(model.load(coefficient_file, parameter_file, elements, error)) << error; + EXPECT_TRUE(model.parameters.switchinner_flag); + ASSERT_EQ(model.parameters.sinner.size(), static_cast(1)); + EXPECT_DOUBLE_EQ(model.parameters.sinner[0], 1.5); + { + std::ofstream output(parameter_file.c_str(), std::ios::app); + output << "unsupported_option 1\n"; + } + EXPECT_FALSE(model.load(coefficient_file, parameter_file, elements, error)); + EXPECT_NE(error.find("unsupported_option"), std::string::npos); + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\ntwojmax 0\nswitchinnerflag 1\nsinner 1.5\n"; + } + error.clear(); + EXPECT_FALSE(model.load(coefficient_file, parameter_file, elements, error)); + EXPECT_NE(error.find("requires sinner and dinner"), std::string::npos); + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\ntwojmax 0\nswitchinnerflag 0\n"; + } + error.clear(); + EXPECT_TRUE(model.load(coefficient_file, parameter_file, elements, error)) << error; + std::remove(coefficient_file.c_str()); + std::remove(parameter_file.c_str()); +} + +TEST(SnapModel, RejectsExtraCoefficientFileData) +{ + const std::string coefficient_file = testing::TempDir() + "snap_extra_coeff.txt"; + const std::string parameter_file = testing::TempDir() + "snap_extra_param.txt"; + { + std::ofstream output(parameter_file.c_str()); + output << "rcutfac 1.0\ntwojmax 0\n"; + } + const std::vector contents = { + "1 2 extra\nSi 2.0 1.0\n0.0\n1.0\n", + "1 2\nSi 2.0 1.0 extra\n0.0\n1.0\n", + "1 2\nSi 2.0 1.0\n0.0 extra\n1.0\n", + "1 2\nSi 2.0 1.0\n0.0\n1.0\nextra\n"}; + const std::vector expected = { + "header", "element line", "coefficient value", "after the SNAP coefficient blocks"}; + for (std::size_t i = 0; i < contents.size(); ++i) + { + { + std::ofstream output(coefficient_file.c_str()); + output << contents[i]; + } + ModuleMD::SNAP::ModelData model; + std::string error; + const std::vector elements(1, "Si"); + EXPECT_FALSE(model.load(coefficient_file, parameter_file, elements, error)); + EXPECT_NE(error.find(expected[i]), std::string::npos) << error; + } + std::remove(coefficient_file.c_str()); + std::remove(parameter_file.c_str()); +} + +} // namespace diff --git a/tests/04_FF/102_SNAP_base/INPUT b/tests/04_FF/102_SNAP_base/INPUT new file mode 100644 index 00000000000..59b27db4ddc --- /dev/null +++ b/tests/04_FF/102_SNAP_base/INPUT @@ -0,0 +1,12 @@ +INPUT_PARAMETERS +calculation md +suffix autotest +esolver_type snap +pot_file ../../SIMULATION_MODELS/snap_102.coeff +snap_param_file ../../SIMULATION_MODELS/snap_102.param +cal_force 1 +cal_stress 1 +md_nstep 1 +md_type nve +md_dt 1 +init_vel 0 diff --git a/tests/04_FF/102_SNAP_base/README b/tests/04_FF/102_SNAP_base/README new file mode 100644 index 00000000000..ed605d2add5 --- /dev/null +++ b/tests/04_FF/102_SNAP_base/README @@ -0,0 +1,14 @@ +CPU SNAP regression case with a two-atom Si periodic structure. + +The example uses `../../SIMULATION_MODELS/snap_102.coeff` for the SNAP +coefficients and `../../SIMULATION_MODELS/snap_102.param` for the SNAP +parameters. + +The coefficient and parameter files are deterministic ABACUS-generated test +fixtures. They are intended only to exercise SNAP parsing, energy, force, and +stress evaluation; they are not transferable or physically validated Si +potentials. + +The reference values are generated by ABACUS and checked by the existing +04_FF regression harness. The large cell keeps periodic image interactions +separate while exercising SNAP energy, force, and stress evaluation. diff --git a/tests/04_FF/102_SNAP_base/STRU b/tests/04_FF/102_SNAP_base/STRU new file mode 100644 index 00000000000..8e59e8a0b97 --- /dev/null +++ b/tests/04_FF/102_SNAP_base/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 28.0855 auto auto + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +Si +0 +2 +0.0 0.0 0.0 m 1 1 1 +2.0 0.0 0.0 m 1 1 1 diff --git a/tests/04_FF/102_SNAP_base/result.ref b/tests/04_FF/102_SNAP_base/result.ref new file mode 100644 index 00000000000..9a188a922ac --- /dev/null +++ b/tests/04_FF/102_SNAP_base/result.ref @@ -0,0 +1,5 @@ +etotref -0.62150213 +etotperatomref -0.3107510650 +totalforceref 5.604078 +totalstressref 4.007952 +totaltimeref 35.66 diff --git a/tests/04_FF/103_SNAP_pbc/INPUT b/tests/04_FF/103_SNAP_pbc/INPUT new file mode 100644 index 00000000000..2eee0601b68 --- /dev/null +++ b/tests/04_FF/103_SNAP_pbc/INPUT @@ -0,0 +1,12 @@ +INPUT_PARAMETERS +calculation md +suffix autotest +esolver_type snap +pot_file ../../SIMULATION_MODELS/snap_103.coeff +snap_param_file ../../SIMULATION_MODELS/snap_103.param +cal_force 1 +cal_stress 1 +md_nstep 1 +md_type nve +md_dt 1 +init_vel 0 diff --git a/tests/04_FF/103_SNAP_pbc/README b/tests/04_FF/103_SNAP_pbc/README new file mode 100644 index 00000000000..f6bf91f480f --- /dev/null +++ b/tests/04_FF/103_SNAP_pbc/README @@ -0,0 +1,9 @@ +CPU SNAP periodic-image regression case. + +This uses the same potential and two-atom separation as 102_SNAP_base, +but places the atoms on opposite sides of the periodic boundary. Their +energy, force magnitude, and stress must match the unwrapped case. The +periodic neighbor contributes through the MDCell ghost-force return path. + +Model files: ../../SIMULATION_MODELS/snap_103.coeff and +../../SIMULATION_MODELS/snap_103.param. diff --git a/tests/04_FF/103_SNAP_pbc/STRU b/tests/04_FF/103_SNAP_pbc/STRU new file mode 100644 index 00000000000..2938b4c6e2b --- /dev/null +++ b/tests/04_FF/103_SNAP_pbc/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 28.0855 auto auto + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +Si +0 +2 +1.0 0.0 0.0 m 1 1 1 +19.0 0.0 0.0 m 1 1 1 diff --git a/tests/04_FF/103_SNAP_pbc/result.ref b/tests/04_FF/103_SNAP_pbc/result.ref new file mode 100644 index 00000000000..9a188a922ac --- /dev/null +++ b/tests/04_FF/103_SNAP_pbc/result.ref @@ -0,0 +1,5 @@ +etotref -0.62150213 +etotperatomref -0.3107510650 +totalforceref 5.604078 +totalstressref 4.007952 +totaltimeref 35.66 diff --git a/tests/04_FF/104_SNAP_inner/INPUT b/tests/04_FF/104_SNAP_inner/INPUT new file mode 100644 index 00000000000..57c29805ee8 --- /dev/null +++ b/tests/04_FF/104_SNAP_inner/INPUT @@ -0,0 +1,12 @@ +INPUT_PARAMETERS +calculation md +suffix autotest +esolver_type snap +pot_file ../../SIMULATION_MODELS/snap_104.coeff +snap_param_file ../../SIMULATION_MODELS/snap_104.param +cal_force 1 +cal_stress 1 +md_nstep 1 +md_type nve +md_dt 1 +init_vel 0 diff --git a/tests/04_FF/104_SNAP_inner/README b/tests/04_FF/104_SNAP_inner/README new file mode 100644 index 00000000000..4cdb3ae4c3b --- /dev/null +++ b/tests/04_FF/104_SNAP_inner/README @@ -0,0 +1,13 @@ +CPU SNAP inner-switch regression case. + +The atoms are inside the smooth inner-switch interval (0.8 to 1.2 Angstrom). +This case verifies model-file parsing and the complete MD energy, force, and +stress path with switchinnerflag enabled. The coefficients are a deterministic +software-test fixture, not a validated physical potential. + +The reference was cross-checked with LAMMPS SNAP (30 Mar 2026 development +build, c46bbef): potential energy -0.24612388 eV, x forces -2.7962 and ++2.7962 eV/Angstrom, and Pxx 3999.6081 bar. ABACUS reports Pxx in kbar. + +Model files: ../../SIMULATION_MODELS/snap_104.coeff and +../../SIMULATION_MODELS/snap_104.param. diff --git a/tests/04_FF/104_SNAP_inner/STRU b/tests/04_FF/104_SNAP_inner/STRU new file mode 100644 index 00000000000..8e59e8a0b97 --- /dev/null +++ b/tests/04_FF/104_SNAP_inner/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Si 28.0855 auto auto + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +Si +0 +2 +0.0 0.0 0.0 m 1 1 1 +2.0 0.0 0.0 m 1 1 1 diff --git a/tests/04_FF/104_SNAP_inner/result.ref b/tests/04_FF/104_SNAP_inner/result.ref new file mode 100644 index 00000000000..c5b2fa6a682 --- /dev/null +++ b/tests/04_FF/104_SNAP_inner/result.ref @@ -0,0 +1,5 @@ +etotref -0.24612388 +etotperatomref -0.1230619400 +totalforceref 5.592410 +totalstressref 3.999606 +totaltimeref 3.91 diff --git a/tests/04_FF/105_SNAP_quad/INPUT b/tests/04_FF/105_SNAP_quad/INPUT new file mode 100644 index 00000000000..6b735a80ab4 --- /dev/null +++ b/tests/04_FF/105_SNAP_quad/INPUT @@ -0,0 +1,12 @@ +INPUT_PARAMETERS +calculation md +suffix autotest +esolver_type snap +pot_file ../../SIMULATION_MODELS/snap_105.coeff +snap_param_file ../../SIMULATION_MODELS/snap_105.param +cal_force 1 +cal_stress 1 +md_nstep 1 +md_type nve +md_dt 1 +init_vel 0 diff --git a/tests/04_FF/105_SNAP_quad/README b/tests/04_FF/105_SNAP_quad/README new file mode 100644 index 00000000000..add668c05ae --- /dev/null +++ b/tests/04_FF/105_SNAP_quad/README @@ -0,0 +1,13 @@ +CPU SNAP regression with three non-collinear atoms, twojmax=2, nonzero +angular coefficients, and nonzero quadratic coefficients. The synthetic +potential is for software testing only. Energy, force, and stress are +cross-checked against LAMMPS SNAP using the same cell and coordinates. + +LAMMPS (30 Mar 2026 development build, c46bbef) gives total potential energy +54.909561 eV, forces (eV/Angstrom) of (67.292804, 57.547644, 0), +(-92.139945, 60.594450, 0), and (24.847141, -118.142093, 0), and Pxx +-122909.08 bar. ABACUS reports stress in kbar. The calculator unit test +checks each force component and compares virial with an energy strain derivative. + +Model files: ../../SIMULATION_MODELS/snap_105.coeff and +../../SIMULATION_MODELS/snap_105.param. diff --git a/tests/04_FF/105_SNAP_quad/STRU b/tests/04_FF/105_SNAP_quad/STRU new file mode 100644 index 00000000000..8a6eaca3fb4 --- /dev/null +++ b/tests/04_FF/105_SNAP_quad/STRU @@ -0,0 +1,20 @@ +ATOMIC_SPECIES +Si 28.0855 auto auto + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +Si +0 +3 +0.0 0.0 0.0 m 1 1 1 +2.0 0.0 0.0 m 1 1 1 +0.5 2.5 0.0 m 1 1 1 diff --git a/tests/04_FF/105_SNAP_quad/result.ref b/tests/04_FF/105_SNAP_quad/result.ref new file mode 100644 index 00000000000..c7bd38b168d --- /dev/null +++ b/tests/04_FF/105_SNAP_quad/result.ref @@ -0,0 +1,5 @@ +etotref 54.909561 +etotperatomref 18.3031870000 +totalforceref 420.564050 +totalstressref 422.994068 +totaltimeref 0.02 diff --git a/tests/04_FF/106_SNAP_scf/INPUT b/tests/04_FF/106_SNAP_scf/INPUT new file mode 100644 index 00000000000..cf8103ce55e --- /dev/null +++ b/tests/04_FF/106_SNAP_scf/INPUT @@ -0,0 +1,12 @@ +INPUT_PARAMETERS +calculation scf +suffix autotest +esolver_type snap +pot_file ../../SIMULATION_MODELS/snap_106.coeff +snap_param_file ../../SIMULATION_MODELS/snap_106.param +cal_force 1 +cal_stress 1 +md_nstep 1 +md_type nve +md_dt 1 +init_vel 0 diff --git a/tests/04_FF/106_SNAP_scf/README b/tests/04_FF/106_SNAP_scf/README new file mode 100644 index 00000000000..932d66466e6 --- /dev/null +++ b/tests/04_FF/106_SNAP_scf/README @@ -0,0 +1,7 @@ +UnitCell SNAP single-point regression using the same angular and quadratic +model as case 105. The calculation enters the UnitCell solver path rather +than the direct MDCell path. Energy, total force, and stress must match case +105 and the LAMMPS reference described there. + +Model files: ../../SIMULATION_MODELS/snap_106.coeff and +../../SIMULATION_MODELS/snap_106.param. diff --git a/tests/04_FF/106_SNAP_scf/STRU b/tests/04_FF/106_SNAP_scf/STRU new file mode 100644 index 00000000000..8a6eaca3fb4 --- /dev/null +++ b/tests/04_FF/106_SNAP_scf/STRU @@ -0,0 +1,20 @@ +ATOMIC_SPECIES +Si 28.0855 auto auto + +LATTICE_CONSTANT +1.0 + +LATTICE_VECTORS +20.0 0.0 0.0 +0.0 20.0 0.0 +0.0 0.0 20.0 + +ATOMIC_POSITIONS +Cartesian + +Si +0 +3 +0.0 0.0 0.0 m 1 1 1 +2.0 0.0 0.0 m 1 1 1 +0.5 2.5 0.0 m 1 1 1 diff --git a/tests/04_FF/106_SNAP_scf/result.ref b/tests/04_FF/106_SNAP_scf/result.ref new file mode 100644 index 00000000000..c7bd38b168d --- /dev/null +++ b/tests/04_FF/106_SNAP_scf/result.ref @@ -0,0 +1,5 @@ +etotref 54.909561 +etotperatomref 18.3031870000 +totalforceref 420.564050 +totalstressref 422.994068 +totaltimeref 0.02 diff --git a/tests/04_FF/CASES_CPU.txt b/tests/04_FF/CASES_CPU.txt index 886fac4f083..4f66a60ae0a 100644 --- a/tests/04_FF/CASES_CPU.txt +++ b/tests/04_FF/CASES_CPU.txt @@ -20,3 +20,8 @@ 20_LJ_dry_run #50_DP_Al #101_NEP_HfO2 +102_SNAP_base +103_SNAP_pbc +104_SNAP_inner +105_SNAP_quad +106_SNAP_scf diff --git a/tests/04_FF/CMakeLists.txt b/tests/04_FF/CMakeLists.txt index 3041d6d389e..b6fed1438da 100644 --- a/tests/04_FF/CMakeLists.txt +++ b/tests/04_FF/CMakeLists.txt @@ -12,3 +12,15 @@ else() COMMAND ${BASH} ../integrate/Autotest.sh -a ${ABACUS_BIN_PATH} -n 4 -j auto WORKING_DIRECTORY ${ABACUS_TEST_DIR}/04_FF) endif() + +find_program(SNAP_TEST_PYTHON python3) +if(ENABLE_MPI AND SNAP_TEST_PYTHON AND MPIEXEC_EXECUTABLE) + add_test( + NAME 04_FF_SNAP_rank_error + COMMAND ${SNAP_TEST_PYTHON} ${CMAKE_CURRENT_SOURCE_DIR}/check_snap_rank_error.py + --abacus ${ABACUS_BIN_PATH} + --mpiexec ${MPIEXEC_EXECUTABLE} + --numproc-flag=${MPIEXEC_NUMPROC_FLAG} + --models ${ABACUS_TEST_DIR}/SIMULATION_MODELS) + set_tests_properties(04_FF_SNAP_rank_error PROPERTIES TIMEOUT 120) +endif() diff --git a/tests/04_FF/check_snap_rank_error.py b/tests/04_FF/check_snap_rank_error.py new file mode 100644 index 00000000000..91ae1bedc6a --- /dev/null +++ b/tests/04_FF/check_snap_rank_error.py @@ -0,0 +1,77 @@ +"""Check that a SNAP error on a nonzero MPI rank remains visible.""" + +import argparse +import os +from pathlib import Path +import subprocess +import tempfile + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--abacus", type=Path, required=True) + parser.add_argument("--mpiexec", type=Path, required=True) + parser.add_argument("--numproc-flag", default="-n") + parser.add_argument("--models", type=Path, required=True) + args = parser.parse_args() + + coefficient = (args.models / "snap_102.coeff").resolve() + with tempfile.TemporaryDirectory(prefix="abacus-snap-rank-") as scratch: + case = Path(scratch) + (case / "INPUT").write_text( + "INPUT_PARAMETERS\n" + "calculation md\n" + "suffix rank_error\n" + "esolver_type snap\n" + f"pot_file {coefficient}\n" + "snap_param_file snap.param\n" + "cal_force 1\n" + "md_nstep 1\n" + "md_type nve\n" + "init_vel 0\n", + encoding="utf-8", + ) + (case / "snap.param").write_text( + "rcutfac 4.9\n" + "twojmax 8\n" + "rfac0 0.99363\n" + "rmin0 0.2\n" + "bzeroflag 0\n", + encoding="utf-8", + ) + # Rank 0 owns the distant atom; rank 1 owns the close pair. + (case / "STRU").write_text( + "ATOMIC_SPECIES\nSi 28.0855 auto auto\n\n" + "LATTICE_CONSTANT\n1.0\n\n" + "LATTICE_VECTORS\n20.0 0.0 0.0\n0.0 20.0 0.0\n0.0 0.0 20.0\n\n" + "ATOMIC_POSITIONS\nCartesian\n\nSi\n0\n3\n" + "2.0 0.0 0.0 m 1 1 1\n" + "14.0 0.0 0.0 m 1 1 1\n" + "14.2 0.0 0.0 m 1 1 1\n", + encoding="utf-8", + ) + environment = os.environ.copy() + environment["OMP_NUM_THREADS"] = "1" + environment["OMPI_ALLOW_RUN_AS_ROOT"] = "1" + environment["OMPI_ALLOW_RUN_AS_ROOT_CONFIRM"] = "1" + command = [str(args.mpiexec), args.numproc_flag, "2", str(args.abacus)] + result = subprocess.run( + command, cwd=str(case), env=environment, text=True, + stdout=subprocess.PIPE, stderr=subprocess.PIPE, timeout=90, + ) + expected = "ESolver_SNAP rank 1: SNAP evaluation failed: " \ + "SNAP neighbor distance must be greater than rmin0." + coefficient_messages = result.stdout.count(" Read SNAP coefficient file: ") + parameter_messages = result.stdout.count(" Read SNAP parameter file: ") + if (result.returncode == 0 or expected not in result.stderr + or "SNAP evaluation failed" in result.stdout + or coefficient_messages != 1 or parameter_messages != 1): + raise AssertionError( + "Expected one rank-0 file confirmation, a rank-1 SNAP error, " + f"and nonzero exit; got {result.returncode}.\n" + f"stdout:\n{result.stdout}\nstderr:\n{result.stderr}" + ) + + +if __name__ == "__main__": + main() diff --git a/tests/SIMULATION_MODELS/snap_102.coeff b/tests/SIMULATION_MODELS/snap_102.coeff new file mode 100644 index 00000000000..fc13f2f7b43 --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_102.coeff @@ -0,0 +1,62 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic coefficient set is created for software testing only. +# It is not a transferable or physically validated material potential. + +1 56 +Si 0.5 1 +0.125000000000 +-0.002000000000 +0.003000000000 +-0.004000000000 +0.005000000000 +-0.006000000000 +0.007000000000 +-0.008000000000 +0.009000000000 +-0.010000000000 +0.011000000000 +-0.012000000000 +0.013000000000 +-0.014000000000 +0.015000000000 +-0.016000000000 +0.017000000000 +-0.018000000000 +0.019000000000 +-0.020000000000 +0.021000000000 +-0.022000000000 +0.023000000000 +-0.024000000000 +0.025000000000 +-0.026000000000 +0.027000000000 +-0.028000000000 +0.029000000000 +-0.030000000000 +0.031000000000 +-0.032000000000 +0.033000000000 +-0.034000000000 +0.035000000000 +-0.036000000000 +0.037000000000 +-0.038000000000 +0.039000000000 +-0.040000000000 +0.041000000000 +-0.042000000000 +0.043000000000 +-0.044000000000 +0.045000000000 +-0.046000000000 +0.047000000000 +-0.048000000000 +0.049000000000 +-0.050000000000 +0.051000000000 +-0.052000000000 +0.053000000000 +-0.054000000000 +0.055000000000 +-0.056000000000 diff --git a/tests/SIMULATION_MODELS/snap_102.param b/tests/SIMULATION_MODELS/snap_102.param new file mode 100644 index 00000000000..23e0d83d38b --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_102.param @@ -0,0 +1,10 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic parameter set is created for software testing only. +# It is not a transferable or physically validated material potential. + +rcutfac 4.9 +twojmax 8 +rfac0 0.99363 +rmin0 0 +quadraticflag 0 +bzeroflag 0 diff --git a/tests/SIMULATION_MODELS/snap_103.coeff b/tests/SIMULATION_MODELS/snap_103.coeff new file mode 100644 index 00000000000..fc13f2f7b43 --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_103.coeff @@ -0,0 +1,62 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic coefficient set is created for software testing only. +# It is not a transferable or physically validated material potential. + +1 56 +Si 0.5 1 +0.125000000000 +-0.002000000000 +0.003000000000 +-0.004000000000 +0.005000000000 +-0.006000000000 +0.007000000000 +-0.008000000000 +0.009000000000 +-0.010000000000 +0.011000000000 +-0.012000000000 +0.013000000000 +-0.014000000000 +0.015000000000 +-0.016000000000 +0.017000000000 +-0.018000000000 +0.019000000000 +-0.020000000000 +0.021000000000 +-0.022000000000 +0.023000000000 +-0.024000000000 +0.025000000000 +-0.026000000000 +0.027000000000 +-0.028000000000 +0.029000000000 +-0.030000000000 +0.031000000000 +-0.032000000000 +0.033000000000 +-0.034000000000 +0.035000000000 +-0.036000000000 +0.037000000000 +-0.038000000000 +0.039000000000 +-0.040000000000 +0.041000000000 +-0.042000000000 +0.043000000000 +-0.044000000000 +0.045000000000 +-0.046000000000 +0.047000000000 +-0.048000000000 +0.049000000000 +-0.050000000000 +0.051000000000 +-0.052000000000 +0.053000000000 +-0.054000000000 +0.055000000000 +-0.056000000000 diff --git a/tests/SIMULATION_MODELS/snap_103.param b/tests/SIMULATION_MODELS/snap_103.param new file mode 100644 index 00000000000..23e0d83d38b --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_103.param @@ -0,0 +1,10 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic parameter set is created for software testing only. +# It is not a transferable or physically validated material potential. + +rcutfac 4.9 +twojmax 8 +rfac0 0.99363 +rmin0 0 +quadraticflag 0 +bzeroflag 0 diff --git a/tests/SIMULATION_MODELS/snap_104.coeff b/tests/SIMULATION_MODELS/snap_104.coeff new file mode 100644 index 00000000000..fc13f2f7b43 --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_104.coeff @@ -0,0 +1,62 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic coefficient set is created for software testing only. +# It is not a transferable or physically validated material potential. + +1 56 +Si 0.5 1 +0.125000000000 +-0.002000000000 +0.003000000000 +-0.004000000000 +0.005000000000 +-0.006000000000 +0.007000000000 +-0.008000000000 +0.009000000000 +-0.010000000000 +0.011000000000 +-0.012000000000 +0.013000000000 +-0.014000000000 +0.015000000000 +-0.016000000000 +0.017000000000 +-0.018000000000 +0.019000000000 +-0.020000000000 +0.021000000000 +-0.022000000000 +0.023000000000 +-0.024000000000 +0.025000000000 +-0.026000000000 +0.027000000000 +-0.028000000000 +0.029000000000 +-0.030000000000 +0.031000000000 +-0.032000000000 +0.033000000000 +-0.034000000000 +0.035000000000 +-0.036000000000 +0.037000000000 +-0.038000000000 +0.039000000000 +-0.040000000000 +0.041000000000 +-0.042000000000 +0.043000000000 +-0.044000000000 +0.045000000000 +-0.046000000000 +0.047000000000 +-0.048000000000 +0.049000000000 +-0.050000000000 +0.051000000000 +-0.052000000000 +0.053000000000 +-0.054000000000 +0.055000000000 +-0.056000000000 diff --git a/tests/SIMULATION_MODELS/snap_104.param b/tests/SIMULATION_MODELS/snap_104.param new file mode 100644 index 00000000000..b73ae7c173b --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_104.param @@ -0,0 +1,13 @@ +# ABACUS-generated SNAP regression fixture +# This deterministic parameter set is created for software testing only. +# It is not a transferable or physically validated material potential. + +rcutfac 4.9 +twojmax 8 +rfac0 0.99363 +rmin0 0 +switchinnerflag 1 +sinner 1.0 +dinner 0.2 +quadraticflag 0 +bzeroflag 0 diff --git a/tests/SIMULATION_MODELS/snap_105.coeff b/tests/SIMULATION_MODELS/snap_105.coeff new file mode 100644 index 00000000000..aa7b790cd3f --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_105.coeff @@ -0,0 +1,24 @@ +# Deterministic software-test coefficients for angular and quadratic SNAP. +1 21 +Si 2.25 1.0 +0.2 +0.3 +-0.1 +0.5 +0.2 +-0.4 +0.25 +0.0 +0.0 +0.0 +-0.15 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 diff --git a/tests/SIMULATION_MODELS/snap_105.param b/tests/SIMULATION_MODELS/snap_105.param new file mode 100644 index 00000000000..eebdca8cf58 --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_105.param @@ -0,0 +1,6 @@ +rcutfac 1.0 +twojmax 2 +rfac0 0.99363 +rmin0 0.0 +quadraticflag 1 +bzeroflag 1 diff --git a/tests/SIMULATION_MODELS/snap_106.coeff b/tests/SIMULATION_MODELS/snap_106.coeff new file mode 100644 index 00000000000..aa7b790cd3f --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_106.coeff @@ -0,0 +1,24 @@ +# Deterministic software-test coefficients for angular and quadratic SNAP. +1 21 +Si 2.25 1.0 +0.2 +0.3 +-0.1 +0.5 +0.2 +-0.4 +0.25 +0.0 +0.0 +0.0 +-0.15 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 +0.0 diff --git a/tests/SIMULATION_MODELS/snap_106.param b/tests/SIMULATION_MODELS/snap_106.param new file mode 100644 index 00000000000..eebdca8cf58 --- /dev/null +++ b/tests/SIMULATION_MODELS/snap_106.param @@ -0,0 +1,6 @@ +rcutfac 1.0 +twojmax 2 +rfac0 0.99363 +rmin0 0.0 +quadraticflag 1 +bzeroflag 1