From 5225369481289bef488061a65c6eed395304ee79 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Tue, 7 Jul 2026 23:26:33 +0800 Subject: [PATCH 01/24] Add fixed-cell i-PI socket driver --- docs/advanced/input_files/input-main.md | 8 +- docs/parameters.yaml | 8 +- source/CMakeLists.txt | 4 +- .../read_input_item_system.cpp | 15 +- .../test_serial/read_input_item_test.cpp | 5 + source/source_main/driver.h | 1 + source/source_main/driver_ipi.cpp | 507 ++++++++++++++++++ source/source_main/driver_run.cpp | 9 +- source/source_main/ipi_socket.cpp | 241 +++++++++ source/source_main/ipi_socket.h | 40 ++ 10 files changed, 829 insertions(+), 9 deletions(-) create mode 100644 source/source_main/driver_ipi.cpp create mode 100644 source/source_main/ipi_socket.cpp create mode 100644 source/source_main/ipi_socket.h diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 334144f44c5..c6c1b732991 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -584,6 +584,7 @@ - relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations - cell-relax: perform cell relaxation calculations - md: perform molecular dynamics simulations + - socket: run as a socket client for external drivers using the i-PI protocol - get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information - get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -840,7 +841,12 @@ ### chg_extrap - **Type**: String -- **Description**: Charge extrapolation method for MD and relaxation calculations. +- **Description**: Charge extrapolation method for MD, relaxation, and socket-driven calculations. + + When set to default, ABACUS chooses second-order for md, first-order for + relax/cell-relax/socket, and atomic for other calculations. Socket-driven + molecular dynamics can explicitly set second-order if the external driver + updates structures smoothly enough for second-order extrapolation. - **Default**: default ### nb2d diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 00ef9b9dcf7..d183783ad81 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -28,6 +28,7 @@ parameters: * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations * cell-relax: perform cell relaxation calculations * md: perform molecular dynamics simulations + * socket: run as a socket client for external drivers using the i-PI protocol * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -328,7 +329,12 @@ parameters: category: System variables type: String description: | - Charge extrapolation method for MD and relaxation calculations. + Charge extrapolation method for MD, relaxation, and socket-driven calculations. + + When set to default, ABACUS chooses second-order for md, first-order for + relax/cell-relax/socket, and atomic for other calculations. Socket-driven + molecular dynamics can explicitly set second-order if the external driver + updates structures smoothly enough for second-order extrapolation. default_value: default unit: "" availability: "" diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 0337707b61a..a5b78bd3b50 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -477,7 +477,9 @@ add_library( driver OBJECT source_main/driver.cpp - source_main/driver_run.cpp) + source_main/driver_run.cpp + source_main/driver_ipi.cpp + source_main/ipi_socket.cpp) list(APPEND device_srcs source_pw/module_pwdft/kernels/nonlocal_op.cpp diff --git a/source/source_io/module_parameter/read_input_item_system.cpp b/source/source_io/module_parameter/read_input_item_system.cpp index 3b0ae011ec4..6ed572b21f9 100644 --- a/source/source_io/module_parameter/read_input_item_system.cpp +++ b/source/source_io/module_parameter/read_input_item_system.cpp @@ -75,7 +75,7 @@ void ReadInput::item_system() } { Input_Item item("calculation"); - item.annotation = "scf; relax; md; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; + item.annotation = "scf; relax; md; socket; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; item.category = "System variables"; item.type = "String"; item.description = R"(Specify the type of calculation. @@ -85,6 +85,7 @@ void ReadInput::item_system() * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations * cell-relax: perform cell relaxation calculations * md: perform molecular dynamics simulations +* socket: run as a socket client for external drivers using the i-PI protocol * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -102,6 +103,7 @@ void ReadInput::item_system() std::vector callist = {"scf", "relax", "md", + "socket", "cell-relax", "nscf", "get_s", @@ -266,7 +268,7 @@ void ReadInput::item_system() item.description = "If set to True, calculate the force at the end of the electronic iteration."; item.default_value = "False"; item.reset_value = [](const Input_Item& item, Parameter& para) { - std::vector use_force = {"cell-relax", "relax", "md"}; + std::vector use_force = {"cell-relax", "relax", "md", "socket"}; std::vector not_use_force = {"get_wf", "get_pchg", "get_s"}; if (std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) { @@ -873,7 +875,12 @@ Available options are: item.annotation = "atomic; first-order; second-order; dm:coefficients of SIA"; item.category = "System variables"; item.type = "String"; - item.description = "Charge extrapolation method for MD and relaxation calculations."; + item.description = R"(Charge extrapolation method for MD, relaxation, and socket-driven calculations. + +When set to default, ABACUS chooses second-order for md, first-order for +relax/cell-relax/socket, and atomic for other calculations. Socket-driven +molecular dynamics can explicitly set second-order if the external driver +updates structures smoothly enough for second-order extrapolation.)"; item.default_value = "default"; read_sync_string(input.chg_extrap); item.reset_value = [](const Input_Item& item, Parameter& para) { @@ -882,7 +889,7 @@ Available options are: para.input.chg_extrap = "second-order"; } else if (para.input.chg_extrap == "default" - && (para.input.calculation == "relax" || para.input.calculation == "cell-relax")) + && (para.input.calculation == "relax" || para.input.calculation == "cell-relax" || para.input.calculation == "socket")) { para.input.chg_extrap = "first-order"; } diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 61673a418c7..eb043fe1b03 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -362,6 +362,11 @@ TEST_F(InputTest, Item_test) it->second.reset_value(it->second, param); EXPECT_EQ(param.input.chg_extrap, "first-order"); + param.input.chg_extrap = "default"; + param.input.calculation = "socket"; + it->second.reset_value(it->second, param); + EXPECT_EQ(param.input.chg_extrap, "first-order"); + param.input.chg_extrap = "default"; param.input.calculation = "none"; it->second.reset_value(it->second, param); diff --git a/source/source_main/driver.h b/source/source_main/driver.h index e0204d079ca..5a51701e26c 100644 --- a/source/source_main/driver.h +++ b/source/source_main/driver.h @@ -36,6 +36,7 @@ class Driver // the actual calculations void driver_run(); + void driver_ipi_run(); // Init harewares according to Input parameters void init_hardware(); diff --git a/source/source_main/driver_ipi.cpp b/source/source_main/driver_ipi.cpp new file mode 100644 index 00000000000..7be2105a9b3 --- /dev/null +++ b/source/source_main/driver_ipi.cpp @@ -0,0 +1,507 @@ +#include "source_main/driver.h" + +#include "source_main/ipi_socket.h" +#include "source_base/global_function.h" +#include "source_base/mathzone.h" +#include "source_cell/check_atomic_stru.h" +#include "source_cell/update_cell.h" +#include "source_esolver/esolver.h" +#include "source_base/global_variable.h" +#include "source_io/module_json/para_json.h" +#include "source_io/module_output/print_info.h" +#include "source_io/module_parameter/parameter.h" + +#ifdef __MPI +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +namespace +{ +constexpr double RY_TO_HARTREE = 0.5; +constexpr int IPI_RANK_ROOT = 0; + +bool is_root() +{ + return GlobalV::MY_RANK == IPI_RANK_ROOT; +} + +void bcast_int(int& value) +{ +#ifdef __MPI + MPI_Bcast(&value, 1, MPI_INT, IPI_RANK_ROOT, MPI_COMM_WORLD); +#else + (void)value; +#endif +} + +void bcast_double_vector(std::vector& values) +{ +#ifdef __MPI + MPI_Bcast(values.data(), static_cast(values.size()), MPI_DOUBLE, IPI_RANK_ROOT, MPI_COMM_WORLD); +#else + (void)values; +#endif +} + +std::string bcast_string(std::string value) +{ + int nbytes = static_cast(value.size()); + bcast_int(nbytes); + if (nbytes < 0) + { + throw std::runtime_error("negative string length in i-PI broadcast"); + } + if (!is_root()) + { + value.assign(static_cast(nbytes), '\0'); + } +#ifdef __MPI + if (nbytes > 0) + { + MPI_Bcast(&value[0], nbytes, MPI_CHAR, IPI_RANK_ROOT, MPI_COMM_WORLD); + } +#endif + return value; +} + +void throw_if_root_io_failed(int root_failed, const std::string& root_message) +{ + bcast_int(root_failed); + const std::string message = bcast_string(root_message); + if (root_failed != 0) + { + throw std::runtime_error(message.empty() ? "i-PI socket I/O failed" : message); + } +} + +std::string bcast_header(const std::string& root_header) +{ + char buffer[13] = {' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', '\0'}; + if (is_root()) + { + const std::size_t n = root_header.size() > 12 ? 12 : root_header.size(); + for (std::size_t i = 0; i < n; ++i) + { + buffer[i] = root_header[i]; + } + } +#ifdef __MPI + MPI_Bcast(buffer, 12, MPI_CHAR, IPI_RANK_ROOT, MPI_COMM_WORLD); +#endif + std::string out(buffer, 12); + while (!out.empty() && out.back() == ' ') + { + out.pop_back(); + } + return out; +} + +std::string ipi_address() +{ + const char* env = std::getenv("ABACUS_IPI_ADDRESS"); + if (env == nullptr || std::string(env).empty()) + { + return "localhost:31415"; + } + return std::string(env); +} + +std::vector ipi_cell_bohr_from_unitcell(const UnitCell& ucell) +{ + const double lat0 = ucell.lat0; + // ASE/i-PI sends POSDATA cell as cell.T in C order. ABACUS stores + // lattice vectors as rows in latvec, so use the transposed order here. + return { + ucell.latvec.e11 * lat0, ucell.latvec.e21 * lat0, ucell.latvec.e31 * lat0, + ucell.latvec.e12 * lat0, ucell.latvec.e22 * lat0, ucell.latvec.e32 * lat0, + ucell.latvec.e13 * lat0, ucell.latvec.e23 * lat0, ucell.latvec.e33 * lat0, + }; +} + +double max_abs_delta(const std::vector& a, const std::vector& b) +{ + if (a.size() != b.size()) + { + return 1.0e99; + } + double out = 0.0; + for (std::size_t i = 0; i < a.size(); ++i) + { + out = std::max(out, std::abs(a[i] - b[i])); + } + return out; +} + +void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& positions_bohr) +{ + if (positions_bohr.size() != static_cast(3 * ucell.nat)) + { + throw std::runtime_error("POSDATA atom count does not match STRU."); + } + int iat = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + Atom* atom = &ucell.atoms[it]; + for (int ia = 0; ia < atom->na; ++ia) + { + const double tau_x = positions_bohr[3 * iat + 0] / ucell.lat0; + const double tau_y = positions_bohr[3 * iat + 1] / ucell.lat0; + const double tau_z = positions_bohr[3 * iat + 2] / ucell.lat0; + + double dx = 0.0; + double dy = 0.0; + double dz = 0.0; + ModuleBase::Mathzone::Cartesian_to_Direct(tau_x, + tau_y, + tau_z, + ucell.latvec.e11, + ucell.latvec.e12, + ucell.latvec.e13, + ucell.latvec.e21, + ucell.latvec.e22, + ucell.latvec.e23, + ucell.latvec.e31, + ucell.latvec.e32, + ucell.latvec.e33, + dx, + dy, + dz); + + atom->dis[ia].x = dx - atom->taud[ia].x; + atom->dis[ia].y = dy - atom->taud[ia].y; + atom->dis[ia].z = dz - atom->taud[ia].z; + atom->taud[ia].x = dx; + atom->taud[ia].y = dy; + atom->taud[ia].z = dz; + atom->tau[ia].x = tau_x; + atom->tau[ia].y = tau_y; + atom->tau[ia].z = tau_z; + ++iat; + } + } + unitcell::periodic_boundary_adjustment(ucell.atoms, ucell.latvec, ucell.ntype); + ucell.ionic_position_updated = true; + ucell.cell_parameter_updated = false; +} + +std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force) +{ + std::vector out(static_cast(force.nr * force.nc)); + for (int iat = 0; iat < force.nr; ++iat) + { + for (int idir = 0; idir < force.nc; ++idir) + { + out[static_cast(3 * iat + idir)] = force(iat, idir) * RY_TO_HARTREE; + } + } + return out; +} + +class CalculationModeGuard +{ + public: + explicit CalculationModeGuard(const std::string& inner_calculation) + : outer_calculation_(PARAM.inp.calculation) + { + const_cast(PARAM.inp.calculation) = inner_calculation; + } + + ~CalculationModeGuard() + { + const_cast(PARAM.inp.calculation) = outer_calculation_; + } + + CalculationModeGuard(const CalculationModeGuard&) = delete; + CalculationModeGuard& operator=(const CalculationModeGuard&) = delete; + + private: + std::string outer_calculation_; +}; +} // namespace + +void Driver::driver_ipi_run() +{ + ModuleBase::TITLE("Driver", "driver_ipi_run"); + + // "socket" is an outer driver mode. The KS/LCAO ESolver internals use the + // standard SCF code path for each POSDATA request from the i-PI protocol. + CalculationModeGuard calculation_guard("scf"); + + UnitCell ucell; + ucell.setup(PARAM.inp.latname, PARAM.inp.ntype, PARAM.inp.lmaxmax, PARAM.inp.init_vel, PARAM.inp.fixed_axes); + ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); + unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); + + IpiSocket socket; + std::unique_ptr p_esolver; + bool hardware_initialized = false; + bool esolver_ready = false; + bool runner_completed = false; + std::string pending_error; + + try + { + this->init_hardware(); + hardware_initialized = true; + + p_esolver.reset(ModuleESolver::init_esolver(PARAM.inp, ucell)); + p_esolver->before_all_runners(ucell, PARAM.inp); + esolver_ready = true; + +#ifdef __RAPIDJSON + Json::gen_stru_wrapper(&ucell); +#endif + + int io_failed = 0; + std::string io_message; + if (is_root()) + { + try + { + const std::string address = ipi_address(); + GlobalV::ofs_running << " ABACUS socket driver connecting to i-PI endpoint " << address << std::endl; + socket.connect(address); + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + + bool isinit = false; + bool hasdata = false; + int istep = 0; + int nat_return = ucell.nat; + double energy_hartree = 0.0; + std::vector forces_hartree_bohr(static_cast(3 * ucell.nat), 0.0); + std::vector virial_hartree(9, 0.0); + + const std::vector reference_cell = ipi_cell_bohr_from_unitcell(ucell); + + while (true) + { + std::string header; + io_failed = 0; + io_message.clear(); + if (is_root()) + { + try + { + header = socket.read_header(); + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + header = bcast_header(header); + + if (header == "STATUS") + { + io_failed = 0; + io_message.clear(); + if (is_root()) + { + try + { + if (hasdata) + { + socket.write_header("HAVEDATA"); + } + else if (isinit) + { + socket.write_header("READY"); + } + else + { + socket.write_header("NEEDINIT"); + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + } + else if (header == "INIT") + { + int rid = 0; + int nbytes = 0; + std::string params; + io_failed = 0; + io_message.clear(); + if (is_root()) + { + try + { + rid = socket.read_int(); + nbytes = socket.read_int(); + if (nbytes < 0) + { + throw std::runtime_error("negative INIT payload length from i-PI socket"); + } + if (nbytes > 0) + { + params = socket.read_string(static_cast(nbytes)); + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + bcast_int(rid); + bcast_int(nbytes); + if (nbytes > 0 && is_root()) + { + GlobalV::ofs_running << " ABACUS socket INIT params " << params << std::endl; + } + isinit = true; + if (is_root()) + { + GlobalV::ofs_running << " ABACUS socket INIT replica " << rid << std::endl; + } + } + else if (header == "POSDATA") + { + std::vector cell(9, 0.0); + std::vector inv_cell(9, 0.0); + int nat_socket = 0; + std::vector positions; + io_failed = 0; + io_message.clear(); + if (is_root()) + { + try + { + cell = socket.read_doubles(9); + inv_cell = socket.read_doubles(9); + nat_socket = socket.read_int(); + if (nat_socket < 0) + { + throw std::runtime_error("negative POSDATA atom count from i-PI socket"); + } + positions = socket.read_doubles(static_cast(3 * nat_socket)); + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + bcast_double_vector(cell); + bcast_double_vector(inv_cell); + bcast_int(nat_socket); + if (!is_root()) + { + positions.assign(static_cast(3 * nat_socket), 0.0); + } + bcast_double_vector(positions); + + if (nat_socket != ucell.nat) + { + throw std::runtime_error("POSDATA atom count does not match STRU."); + } + const double max_cell_delta_bohr = max_abs_delta(cell, reference_cell); + if (max_cell_delta_bohr > 1.0e-6) + { + throw std::runtime_error("variable-cell socket updates are not supported yet."); + } + + set_positions_from_ipi_bohr(ucell, positions); + p_esolver->runner(ucell, istep); + runner_completed = true; + energy_hartree = p_esolver->cal_energy() * RY_TO_HARTREE; + ModuleBase::matrix force; + if (PARAM.inp.cal_force) + { + p_esolver->cal_force(ucell, force); + forces_hartree_bohr = flatten_forces_hartree_per_bohr(force); + } + ++istep; + hasdata = true; + } + else if (header == "GETFORCE") + { + io_failed = 0; + io_message.clear(); + if (is_root()) + { + try + { + socket.write_header("FORCEREADY"); + socket.write_double(energy_hartree); + socket.write_int(nat_return); + socket.write_doubles(forces_hartree_bohr); + socket.write_doubles(virial_hartree); + socket.write_int(0); + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + throw_if_root_io_failed(io_failed, io_message); + isinit = false; + hasdata = false; + } + else + { + if (is_root()) + { + GlobalV::ofs_running << " ABACUS socket driver exiting on header " << header << std::endl; + } + break; + } + } + } + catch (const std::exception& exc) + { + pending_error = exc.what(); + if (is_root()) + { + GlobalV::ofs_running << " ABACUS socket driver ended with error: " << pending_error << std::endl; + } + } + + if (is_root()) + { + socket.close(); + } + if (esolver_ready && runner_completed && p_esolver) + { + p_esolver->after_all_runners(ucell); + } + p_esolver.reset(); + if (hardware_initialized) + { + this->finalize_hardware(); + } + +#ifdef __RAPIDJSON + Json::create_Json(&ucell, PARAM); +#endif + + if (!pending_error.empty()) + { + ModuleBase::WARNING_QUIT("ABACUS socket", pending_error); + } +} diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 4911b133de0..9285f9b2330 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -38,6 +38,13 @@ void Driver::driver_run() { ModuleBase::TITLE("Driver", "driver_run"); + const std::string cal = PARAM.inp.calculation; + if (cal == "socket") + { + this->driver_ipi_run(); + return; + } + //! 1: setup cell and atom information // this warning should not be here, mohan 2024-05-22 #ifndef __LCAO @@ -71,8 +78,6 @@ void Driver::driver_run() Json::gen_stru_wrapper(&ucell); #endif - const std::string cal = PARAM.inp.calculation; - //! 4: different types of calculations if (cal == "md") { diff --git a/source/source_main/ipi_socket.cpp b/source/source_main/ipi_socket.cpp new file mode 100644 index 00000000000..eca053566ff --- /dev/null +++ b/source/source_main/ipi_socket.cpp @@ -0,0 +1,241 @@ +#include "ipi_socket.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ +constexpr std::size_t IPI_HEADER_LEN = 12; + +std::string errno_message(const std::string& prefix) +{ + return prefix + ": " + std::strerror(errno); +} + +std::string trim_header(const char* data) +{ + std::string value(data, IPI_HEADER_LEN); + while (!value.empty() && value.back() == ' ') + { + value.pop_back(); + } + return value; +} + +std::string padded_header(const std::string& header) +{ + if (header.size() > IPI_HEADER_LEN) + { + throw std::runtime_error("i-PI header is longer than 12 bytes: " + header); + } + std::string out = header; + out.resize(IPI_HEADER_LEN, ' '); + return out; +} +} // namespace + +IpiSocket::~IpiSocket() +{ + this->close(); +} + +void IpiSocket::connect(const std::string& address) +{ + this->close(); + const std::size_t colon = address.rfind(':'); + if (colon == std::string::npos) + { + throw std::runtime_error("i-PI address must be host:port or path:UNIX, got " + address); + } + const std::string host = address.substr(0, colon); + const std::string service = address.substr(colon + 1); + + if (service == "UNIX") + { + fd_ = ::socket(AF_UNIX, SOCK_STREAM, 0); + if (fd_ < 0) + { + throw std::runtime_error(errno_message("failed to create UNIX socket")); + } + sockaddr_un addr; + std::memset(&addr, 0, sizeof(addr)); + addr.sun_family = AF_UNIX; + if (host.size() >= sizeof(addr.sun_path)) + { + this->close(); + throw std::runtime_error("UNIX socket path too long: " + host); + } + std::strncpy(addr.sun_path, host.c_str(), sizeof(addr.sun_path) - 1); + if (::connect(fd_, reinterpret_cast(&addr), sizeof(addr)) != 0) + { + const std::string msg = errno_message("failed to connect UNIX i-PI socket " + host); + this->close(); + throw std::runtime_error(msg); + } + return; + } + + addrinfo hints; + std::memset(&hints, 0, sizeof(hints)); + hints.ai_family = AF_UNSPEC; + hints.ai_socktype = SOCK_STREAM; + + addrinfo* result = nullptr; + const int gai = ::getaddrinfo(host.c_str(), service.c_str(), &hints, &result); + if (gai != 0) + { + throw std::runtime_error("failed to resolve i-PI socket " + address + ": " + ::gai_strerror(gai)); + } + + std::string last_error; + for (addrinfo* rp = result; rp != nullptr; rp = rp->ai_next) + { + fd_ = ::socket(rp->ai_family, rp->ai_socktype, rp->ai_protocol); + if (fd_ < 0) + { + last_error = errno_message("failed to create INET socket"); + continue; + } + if (::connect(fd_, rp->ai_addr, rp->ai_addrlen) == 0) + { + ::freeaddrinfo(result); + return; + } + last_error = errno_message("failed to connect INET i-PI socket " + address); + this->close(); + } + ::freeaddrinfo(result); + throw std::runtime_error(last_error.empty() ? "failed to connect i-PI socket " + address : last_error); +} + +void IpiSocket::close() +{ + if (fd_ >= 0) + { + ::close(fd_); + fd_ = -1; + } +} + +std::string IpiSocket::read_header() +{ + char header[IPI_HEADER_LEN]; + this->read_exact(header, sizeof(header)); + return trim_header(header); +} + +void IpiSocket::write_header(const std::string& header) +{ + const std::string padded = padded_header(header); + this->write_exact(padded.data(), padded.size()); +} + +int IpiSocket::read_int() +{ + int value = 0; + this->read_exact(&value, sizeof(value)); + return value; +} + +void IpiSocket::write_int(int value) +{ + this->write_exact(&value, sizeof(value)); +} + +double IpiSocket::read_double() +{ + double value = 0.0; + this->read_exact(&value, sizeof(value)); + return value; +} + +void IpiSocket::write_double(double value) +{ + this->write_exact(&value, sizeof(value)); +} + +std::vector IpiSocket::read_doubles(std::size_t n) +{ + std::vector values(n); + if (!values.empty()) + { + this->read_exact(values.data(), values.size() * sizeof(double)); + } + return values; +} + +void IpiSocket::write_doubles(const std::vector& values) +{ + if (!values.empty()) + { + this->write_exact(values.data(), values.size() * sizeof(double)); + } +} + +std::string IpiSocket::read_string(std::size_t nbytes) +{ + std::string value(nbytes, '\0'); + if (nbytes > 0) + { + this->read_exact(&value[0], nbytes); + } + return value; +} + +void IpiSocket::read_exact(void* data, std::size_t nbytes) +{ + char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { + const ssize_t nread = ::recv(fd_, cursor + done, nbytes - done, 0); + if (nread == 0) + { + throw std::runtime_error("i-PI socket closed while reading"); + } + if (nread < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("i-PI socket read failed")); + } + done += static_cast(nread); + } +} + +void IpiSocket::write_exact(const void* data, std::size_t nbytes) +{ + const char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { +#ifdef MSG_NOSIGNAL + const int flags = MSG_NOSIGNAL; +#else + const int flags = 0; +#endif + const ssize_t nwritten = ::send(fd_, cursor + done, nbytes - done, flags); + if (nwritten == 0) + { + throw std::runtime_error("i-PI socket closed while writing"); + } + if (nwritten < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("i-PI socket write failed")); + } + done += static_cast(nwritten); + } +} diff --git a/source/source_main/ipi_socket.h b/source/source_main/ipi_socket.h new file mode 100644 index 00000000000..3e1713e868d --- /dev/null +++ b/source/source_main/ipi_socket.h @@ -0,0 +1,40 @@ +#ifndef ABACUS_IPI_SOCKET_H +#define ABACUS_IPI_SOCKET_H + +#include +#include +#include + +class IpiSocket +{ + public: + IpiSocket() = default; + ~IpiSocket(); + + IpiSocket(const IpiSocket&) = delete; + IpiSocket& operator=(const IpiSocket&) = delete; + + void connect(const std::string& address); + void close(); + + std::string read_header(); + void write_header(const std::string& header); + + int read_int(); + void write_int(int value); + + double read_double(); + void write_double(double value); + + std::vector read_doubles(std::size_t n); + void write_doubles(const std::vector& values); + std::string read_string(std::size_t nbytes); + + private: + int fd_ = -1; + + void read_exact(void* data, std::size_t nbytes); + void write_exact(const void* data, std::size_t nbytes); +}; + +#endif From e0f70ca031422c4821d5304050c0035f1be85ec7 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Thu, 9 Jul 2026 09:13:53 +0800 Subject: [PATCH 02/24] Refactor i-PI socket driver into source_relax Route calculation=socket through Relax_Driver, keep ESolver lifecycle in driver_run, and add source_relax socket transport tests. --- source/CMakeLists.txt | 4 +- source/source_main/driver.h | 1 - source/source_main/driver_run.cpp | 22 +- source/source_relax/CMakeLists.txt | 2 + .../ipi_socket.cpp | 31 ++- .../ipi_socket.h | 7 + source/source_relax/relax_driver.cpp | 9 + .../socket_driver.cpp} | 243 ++++++----------- source/source_relax/socket_driver.h | 22 ++ source/source_relax/test/CMakeLists.txt | 6 + source/source_relax/test/ipi_socket_test.cpp | 256 ++++++++++++++++++ 11 files changed, 431 insertions(+), 172 deletions(-) rename source/{source_main => source_relax}/ipi_socket.cpp (87%) rename source/{source_main => source_relax}/ipi_socket.h (85%) rename source/{source_main/driver_ipi.cpp => source_relax/socket_driver.cpp} (64%) create mode 100644 source/source_relax/socket_driver.h create mode 100644 source/source_relax/test/ipi_socket_test.cpp diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index a5b78bd3b50..0337707b61a 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -477,9 +477,7 @@ add_library( driver OBJECT source_main/driver.cpp - source_main/driver_run.cpp - source_main/driver_ipi.cpp - source_main/ipi_socket.cpp) + source_main/driver_run.cpp) list(APPEND device_srcs source_pw/module_pwdft/kernels/nonlocal_op.cpp diff --git a/source/source_main/driver.h b/source/source_main/driver.h index 5a51701e26c..e0204d079ca 100644 --- a/source/source_main/driver.h +++ b/source/source_main/driver.h @@ -36,7 +36,6 @@ class Driver // the actual calculations void driver_run(); - void driver_ipi_run(); // Init harewares according to Input parameters void init_hardware(); diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 9285f9b2330..3e7d669b6c2 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -39,11 +39,7 @@ void Driver::driver_run() ModuleBase::TITLE("Driver", "driver_run"); const std::string cal = PARAM.inp.calculation; - if (cal == "socket") - { - this->driver_ipi_run(); - return; - } + const bool socket_mode = (cal == "socket"); //! 1: setup cell and atom information // this warning should not be here, mohan 2024-05-22 @@ -66,12 +62,20 @@ void Driver::driver_run() unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) + Input_para esolver_inp = PARAM.inp; + Input_para driver_inp = PARAM.inp; + if (socket_mode) + { + esolver_inp.calculation = "scf"; + driver_inp.calculation = cal; + } + this->init_hardware(); - ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); + ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(esolver_inp, ucell); //! 3: initialize Esolver and fill json-structure - p_esolver->before_all_runners(ucell, PARAM.inp); + p_esolver->before_all_runners(ucell, esolver_inp); // this Json part should be moved to before_all_runners, mohan 2024-05-12 #ifdef __RAPIDJSON @@ -83,10 +87,10 @@ void Driver::driver_run() { Run_MD::md_line(ucell, p_esolver, PARAM); } - else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf") + else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf" || cal == "socket") { Relax_Driver rl_driver; - rl_driver.relax_driver(p_esolver, ucell, PARAM.inp, GlobalV::ofs_running); + rl_driver.relax_driver(p_esolver, ucell, driver_inp, GlobalV::ofs_running); } else if (cal == "get_s") { diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 9e7ef96b0e2..84f69c06c80 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -2,6 +2,8 @@ add_library( relax OBJECT relax_data.cpp + ipi_socket.cpp + socket_driver.cpp cg_base.cpp relax_driver.cpp relax_sync.cpp diff --git a/source/source_main/ipi_socket.cpp b/source/source_relax/ipi_socket.cpp similarity index 87% rename from source/source_main/ipi_socket.cpp rename to source/source_relax/ipi_socket.cpp index eca053566ff..304a0181813 100644 --- a/source/source_main/ipi_socket.cpp +++ b/source/source_relax/ipi_socket.cpp @@ -1,4 +1,4 @@ -#include "ipi_socket.h" +#include "source_relax/ipi_socket.h" #include #include @@ -41,6 +41,10 @@ std::string padded_header(const std::string& header) } } // namespace +IpiSocketClosed::IpiSocketClosed(const std::string& message) : std::runtime_error(message) +{ +} + IpiSocket::~IpiSocket() { this->close(); @@ -127,7 +131,28 @@ void IpiSocket::close() std::string IpiSocket::read_header() { char header[IPI_HEADER_LEN]; - this->read_exact(header, sizeof(header)); + std::size_t done = 0; + while (done < sizeof(header)) + { + const ssize_t nread = ::recv(fd_, header + done, sizeof(header) - done, 0); + if (nread == 0) + { + if (done == 0) + { + throw IpiSocketClosed("i-PI socket closed before next header"); + } + throw std::runtime_error("i-PI socket closed while reading header"); + } + if (nread < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("i-PI socket header read failed")); + } + done += static_cast(nread); + } return trim_header(header); } @@ -198,7 +223,7 @@ void IpiSocket::read_exact(void* data, std::size_t nbytes) const ssize_t nread = ::recv(fd_, cursor + done, nbytes - done, 0); if (nread == 0) { - throw std::runtime_error("i-PI socket closed while reading"); + throw IpiSocketClosed("i-PI socket closed while reading"); } if (nread < 0) { diff --git a/source/source_main/ipi_socket.h b/source/source_relax/ipi_socket.h similarity index 85% rename from source/source_main/ipi_socket.h rename to source/source_relax/ipi_socket.h index 3e1713e868d..28ce2cb81ff 100644 --- a/source/source_main/ipi_socket.h +++ b/source/source_relax/ipi_socket.h @@ -2,9 +2,16 @@ #define ABACUS_IPI_SOCKET_H #include +#include #include #include +class IpiSocketClosed : public std::runtime_error +{ + public: + explicit IpiSocketClosed(const std::string& message); +}; + class IpiSocket { public: diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 26128bf0408..0cada1bfc3e 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -1,4 +1,5 @@ #include "relax_driver.h" +#include "socket_driver.h" #include "source_base/global_file.h" #include "source_io/module_output/cif_io.h" #include "source_io/module_json/output_info.h" @@ -17,6 +18,14 @@ void Relax_Driver::relax_driver( ModuleBase::TITLE("Relax_Driver", "relax_driver"); ModuleBase::timer::start("Relax_Driver", "relax_driver"); + if (inp.calculation == "socket") + { + Socket_Driver socket_driver; + socket_driver.socket_driver(p_esolver, ucell, inp, ofs_running); + ModuleBase::timer::end("Relax_Driver", "relax_driver"); + return; + } + this->init_relax(ucell.nat, inp); // steps[0]: istep (main iteration step) diff --git a/source/source_main/driver_ipi.cpp b/source/source_relax/socket_driver.cpp similarity index 64% rename from source/source_main/driver_ipi.cpp rename to source/source_relax/socket_driver.cpp index 7be2105a9b3..02928fec9d3 100644 --- a/source/source_main/driver_ipi.cpp +++ b/source/source_relax/socket_driver.cpp @@ -1,25 +1,17 @@ -#include "source_main/driver.h" +#include "socket_driver.h" -#include "source_main/ipi_socket.h" +#include "source_relax/ipi_socket.h" #include "source_base/global_function.h" +#include "source_base/global_variable.h" #include "source_base/mathzone.h" -#include "source_cell/check_atomic_stru.h" +#include "source_base/parallel_common.h" +#include "source_base/timer.h" #include "source_cell/update_cell.h" -#include "source_esolver/esolver.h" -#include "source_base/global_variable.h" -#include "source_io/module_json/para_json.h" -#include "source_io/module_output/print_info.h" -#include "source_io/module_parameter/parameter.h" - -#ifdef __MPI -#include -#endif #include #include #include -#include -#include +#include #include #include @@ -33,75 +25,42 @@ bool is_root() return GlobalV::MY_RANK == IPI_RANK_ROOT; } -void bcast_int(int& value) -{ -#ifdef __MPI - MPI_Bcast(&value, 1, MPI_INT, IPI_RANK_ROOT, MPI_COMM_WORLD); -#else - (void)value; -#endif -} - void bcast_double_vector(std::vector& values) { -#ifdef __MPI - MPI_Bcast(values.data(), static_cast(values.size()), MPI_DOUBLE, IPI_RANK_ROOT, MPI_COMM_WORLD); -#else - (void)values; -#endif + if (!values.empty()) + { + Parallel_Common::bcast_double(values.data(), static_cast(values.size())); + } } -std::string bcast_string(std::string value) +void bcast_socket_string(std::string& value) { - int nbytes = static_cast(value.size()); - bcast_int(nbytes); - if (nbytes < 0) - { - throw std::runtime_error("negative string length in i-PI broadcast"); - } + int size = static_cast(value.size()); + Parallel_Common::bcast_int(size); if (!is_root()) { - value.assign(static_cast(nbytes), '\0'); + value.resize(static_cast(size)); } -#ifdef __MPI - if (nbytes > 0) + if (size > 0) { - MPI_Bcast(&value[0], nbytes, MPI_CHAR, IPI_RANK_ROOT, MPI_COMM_WORLD); + Parallel_Common::bcast_char(&value[0], size); } -#endif - return value; } -void throw_if_root_io_failed(int root_failed, const std::string& root_message) +void quit_if_root_io_failed(int root_failed, std::string root_message) { - bcast_int(root_failed); - const std::string message = bcast_string(root_message); + Parallel_Common::bcast_int(root_failed); + bcast_socket_string(root_message); if (root_failed != 0) { - throw std::runtime_error(message.empty() ? "i-PI socket I/O failed" : message); + ModuleBase::WARNING_QUIT("ABACUS socket", root_message.empty() ? "i-PI socket I/O failed" : root_message); } } -std::string bcast_header(const std::string& root_header) +std::string bcast_header(std::string header) { - char buffer[13] = {' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', ' ', '\0'}; - if (is_root()) - { - const std::size_t n = root_header.size() > 12 ? 12 : root_header.size(); - for (std::size_t i = 0; i < n; ++i) - { - buffer[i] = root_header[i]; - } - } -#ifdef __MPI - MPI_Bcast(buffer, 12, MPI_CHAR, IPI_RANK_ROOT, MPI_COMM_WORLD); -#endif - std::string out(buffer, 12); - while (!out.empty() && out.back() == ' ') - { - out.pop_back(); - } - return out; + bcast_socket_string(header); + return header; } std::string ipi_address() @@ -144,8 +103,9 @@ void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& pos { if (positions_bohr.size() != static_cast(3 * ucell.nat)) { - throw std::runtime_error("POSDATA atom count does not match STRU."); + ModuleBase::WARNING_QUIT("ABACUS socket", "POSDATA atom count does not match STRU."); } + int iat = 0; for (int it = 0; it < ucell.ntype; ++it) { @@ -204,62 +164,29 @@ std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& fo } return out; } +} // namespace -class CalculationModeGuard +void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + std::ofstream& ofs_running) { - public: - explicit CalculationModeGuard(const std::string& inner_calculation) - : outer_calculation_(PARAM.inp.calculation) + ModuleBase::TITLE("Socket_Driver", "socket_driver"); + ModuleBase::timer::start("Socket_Driver", "socket_driver"); + + if (p_esolver == nullptr) { - const_cast(PARAM.inp.calculation) = inner_calculation; + ModuleBase::WARNING_QUIT("ABACUS socket", "socket driver requires a valid ESolver."); } - - ~CalculationModeGuard() + if (!inp.cal_force) { - const_cast(PARAM.inp.calculation) = outer_calculation_; + ModuleBase::WARNING_QUIT("ABACUS socket", "socket calculation requires cal_force=1 for i-PI GETFORCE."); } - CalculationModeGuard(const CalculationModeGuard&) = delete; - CalculationModeGuard& operator=(const CalculationModeGuard&) = delete; - - private: - std::string outer_calculation_; -}; -} // namespace - -void Driver::driver_ipi_run() -{ - ModuleBase::TITLE("Driver", "driver_ipi_run"); - - // "socket" is an outer driver mode. The KS/LCAO ESolver internals use the - // standard SCF code path for each POSDATA request from the i-PI protocol. - CalculationModeGuard calculation_guard("scf"); - - UnitCell ucell; - ucell.setup(PARAM.inp.latname, PARAM.inp.ntype, PARAM.inp.lmaxmax, PARAM.inp.init_vel, PARAM.inp.fixed_axes); - ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); - unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); - IpiSocket socket; - std::unique_ptr p_esolver; - bool hardware_initialized = false; - bool esolver_ready = false; - bool runner_completed = false; - std::string pending_error; try { - this->init_hardware(); - hardware_initialized = true; - - p_esolver.reset(ModuleESolver::init_esolver(PARAM.inp, ucell)); - p_esolver->before_all_runners(ucell, PARAM.inp); - esolver_ready = true; - -#ifdef __RAPIDJSON - Json::gen_stru_wrapper(&ucell); -#endif - int io_failed = 0; std::string io_message; if (is_root()) @@ -267,7 +194,7 @@ void Driver::driver_ipi_run() try { const std::string address = ipi_address(); - GlobalV::ofs_running << " ABACUS socket driver connecting to i-PI endpoint " << address << std::endl; + ofs_running << " ABACUS socket driver connecting to i-PI endpoint " << address << std::endl; socket.connect(address); } catch (const std::exception& exc) @@ -276,12 +203,12 @@ void Driver::driver_ipi_run() io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); + quit_if_root_io_failed(io_failed, io_message); bool isinit = false; bool hasdata = false; int istep = 0; - int nat_return = ucell.nat; + const int nat_return = ucell.nat; double energy_hartree = 0.0; std::vector forces_hartree_bohr(static_cast(3 * ucell.nat), 0.0); std::vector virial_hartree(9, 0.0); @@ -299,16 +226,28 @@ void Driver::driver_ipi_run() { header = socket.read_header(); } + catch (const IpiSocketClosed&) + { + header.clear(); + } catch (const std::exception& exc) { io_failed = 1; io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); + quit_if_root_io_failed(io_failed, io_message); header = bcast_header(header); - if (header == "STATUS") + if (header.empty()) + { + if (is_root()) + { + ofs_running << " ABACUS socket driver exiting after peer closed connection" << std::endl; + } + break; + } + else if (header == "STATUS") { io_failed = 0; io_message.clear(); @@ -335,7 +274,7 @@ void Driver::driver_ipi_run() io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); + quit_if_root_io_failed(io_failed, io_message); } else if (header == "INIT") { @@ -352,9 +291,10 @@ void Driver::driver_ipi_run() nbytes = socket.read_int(); if (nbytes < 0) { - throw std::runtime_error("negative INIT payload length from i-PI socket"); + io_failed = 1; + io_message = "negative INIT payload length from i-PI socket"; } - if (nbytes > 0) + else if (nbytes > 0) { params = socket.read_string(static_cast(nbytes)); } @@ -365,17 +305,17 @@ void Driver::driver_ipi_run() io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); - bcast_int(rid); - bcast_int(nbytes); + quit_if_root_io_failed(io_failed, io_message); + Parallel_Common::bcast_int(rid); + Parallel_Common::bcast_int(nbytes); if (nbytes > 0 && is_root()) { - GlobalV::ofs_running << " ABACUS socket INIT params " << params << std::endl; + ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; } isinit = true; if (is_root()) { - GlobalV::ofs_running << " ABACUS socket INIT replica " << rid << std::endl; + ofs_running << " ABACUS socket INIT replica " << rid << std::endl; } } else if (header == "POSDATA") @@ -395,9 +335,13 @@ void Driver::driver_ipi_run() nat_socket = socket.read_int(); if (nat_socket < 0) { - throw std::runtime_error("negative POSDATA atom count from i-PI socket"); + io_failed = 1; + io_message = "negative POSDATA atom count from i-PI socket"; + } + else + { + positions = socket.read_doubles(static_cast(3 * nat_socket)); } - positions = socket.read_doubles(static_cast(3 * nat_socket)); } catch (const std::exception& exc) { @@ -405,10 +349,10 @@ void Driver::driver_ipi_run() io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); + quit_if_root_io_failed(io_failed, io_message); bcast_double_vector(cell); bcast_double_vector(inv_cell); - bcast_int(nat_socket); + Parallel_Common::bcast_int(nat_socket); if (!is_root()) { positions.assign(static_cast(3 * nat_socket), 0.0); @@ -417,20 +361,27 @@ void Driver::driver_ipi_run() if (nat_socket != ucell.nat) { - throw std::runtime_error("POSDATA atom count does not match STRU."); + ModuleBase::WARNING_QUIT("ABACUS socket", "POSDATA atom count does not match STRU."); } const double max_cell_delta_bohr = max_abs_delta(cell, reference_cell); if (max_cell_delta_bohr > 1.0e-6) { - throw std::runtime_error("variable-cell socket updates are not supported yet."); + ModuleBase::WARNING_QUIT("ABACUS socket", "variable-cell socket updates are not supported yet."); } set_positions_from_ipi_bohr(ucell, positions); p_esolver->runner(ucell, istep); - runner_completed = true; - energy_hartree = p_esolver->cal_energy() * RY_TO_HARTREE; + const double energy_ry = p_esolver->cal_energy(); + energy_hartree = energy_ry * RY_TO_HARTREE; + if (is_root()) + { + ofs_running << " ABACUS socket return energy " + << energy_ry << " Ry, " + << energy_ry * ModuleBase::Ry_to_eV << " eV, " + << energy_hartree << " Ha" << std::endl; + } ModuleBase::matrix force; - if (PARAM.inp.cal_force) + if (inp.cal_force) { p_esolver->cal_force(ucell, force); forces_hartree_bohr = flatten_forces_hartree_per_bohr(force); @@ -459,7 +410,7 @@ void Driver::driver_ipi_run() io_message = exc.what(); } } - throw_if_root_io_failed(io_failed, io_message); + quit_if_root_io_failed(io_failed, io_message); isinit = false; hasdata = false; } @@ -467,7 +418,7 @@ void Driver::driver_ipi_run() { if (is_root()) { - GlobalV::ofs_running << " ABACUS socket driver exiting on header " << header << std::endl; + ofs_running << " ABACUS socket driver exiting on header " << header << std::endl; } break; } @@ -475,33 +426,13 @@ void Driver::driver_ipi_run() } catch (const std::exception& exc) { - pending_error = exc.what(); - if (is_root()) - { - GlobalV::ofs_running << " ABACUS socket driver ended with error: " << pending_error << std::endl; - } + ModuleBase::WARNING_QUIT("ABACUS socket", exc.what()); } if (is_root()) { socket.close(); } - if (esolver_ready && runner_completed && p_esolver) - { - p_esolver->after_all_runners(ucell); - } - p_esolver.reset(); - if (hardware_initialized) - { - this->finalize_hardware(); - } -#ifdef __RAPIDJSON - Json::create_Json(&ucell, PARAM); -#endif - - if (!pending_error.empty()) - { - ModuleBase::WARNING_QUIT("ABACUS socket", pending_error); - } + ModuleBase::timer::end("Socket_Driver", "socket_driver"); } diff --git a/source/source_relax/socket_driver.h b/source/source_relax/socket_driver.h new file mode 100644 index 00000000000..51b75289cb3 --- /dev/null +++ b/source/source_relax/socket_driver.h @@ -0,0 +1,22 @@ +#ifndef ABACUS_SOURCE_RELAX_SOCKET_DRIVER_H +#define ABACUS_SOURCE_RELAX_SOCKET_DRIVER_H + +#include "source_cell/unitcell.h" +#include "source_esolver/esolver.h" +#include "source_io/module_parameter/input_parameter.h" + +#include + +class Socket_Driver +{ + public: + Socket_Driver() = default; + ~Socket_Driver() = default; + + void socket_driver(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + std::ofstream& ofs_running); +}; + +#endif diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 2262fefb419..d1c867dc690 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -6,6 +6,12 @@ abacus_disable_feature_definitions(__ROCM) install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +AddTest( + TARGET MODULE_RELAX_ipi_socket_test + SOURCES ipi_socket_test.cpp ../ipi_socket.cpp +) + AddTest( TARGET MODULE_RELAX_relax_new_line_search LIBS parameter diff --git a/source/source_relax/test/ipi_socket_test.cpp b/source/source_relax/test/ipi_socket_test.cpp new file mode 100644 index 00000000000..21db9be2b90 --- /dev/null +++ b/source/source_relax/test/ipi_socket_test.cpp @@ -0,0 +1,256 @@ +#include "../ipi_socket.h" + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ +constexpr std::size_t IPI_HEADER_LEN = 12; + +std::string errno_message(const std::string& prefix) +{ + return prefix + ": " + std::strerror(errno); +} + +void send_all(int fd, const void* data, std::size_t nbytes) +{ + const char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { + const ssize_t sent = ::send(fd, cursor + done, nbytes - done, 0); + if (sent < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("send failed")); + } + if (sent == 0) + { + throw std::runtime_error("send returned zero"); + } + done += static_cast(sent); + } +} + +void recv_all(int fd, void* data, std::size_t nbytes) +{ + char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { + const ssize_t received = ::recv(fd, cursor + done, nbytes - done, 0); + if (received < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("recv failed")); + } + if (received == 0) + { + throw std::runtime_error("socket closed while receiving test data"); + } + done += static_cast(received); + } +} + +std::string padded_header(const std::string& header) +{ + std::string padded = header; + padded.resize(IPI_HEADER_LEN, ' '); + return padded; +} + +class UnixSocketServer +{ + public: + UnixSocketServer() + { + char dir_template[] = "/tmp/abacus_ipi_socket_test_XXXXXX"; + char* made_dir = ::mkdtemp(dir_template); + if (made_dir == nullptr) + { + throw std::runtime_error(errno_message("mkdtemp failed")); + } + dir_ = made_dir; + path_ = dir_ + "/ipi.sock"; + + listen_fd_ = ::socket(AF_UNIX, SOCK_STREAM, 0); + if (listen_fd_ < 0) + { + throw std::runtime_error(errno_message("socket failed")); + } + + sockaddr_un addr; + std::memset(&addr, 0, sizeof(addr)); + addr.sun_family = AF_UNIX; + std::strncpy(addr.sun_path, path_.c_str(), sizeof(addr.sun_path) - 1); + if (::bind(listen_fd_, reinterpret_cast(&addr), sizeof(addr)) != 0) + { + throw std::runtime_error(errno_message("bind failed")); + } + if (::listen(listen_fd_, 1) != 0) + { + throw std::runtime_error(errno_message("listen failed")); + } + } + + ~UnixSocketServer() + { + if (listen_fd_ >= 0) + { + ::close(listen_fd_); + } + if (!path_.empty()) + { + ::unlink(path_.c_str()); + } + if (!dir_.empty()) + { + ::rmdir(dir_.c_str()); + } + } + + UnixSocketServer(const UnixSocketServer&) = delete; + UnixSocketServer& operator=(const UnixSocketServer&) = delete; + + std::string address() const + { + return path_ + ":UNIX"; + } + + int accept_once() + { + const int fd = ::accept(listen_fd_, nullptr, nullptr); + if (fd < 0) + { + throw std::runtime_error(errno_message("accept failed")); + } + return fd; + } + + private: + int listen_fd_ = -1; + std::string dir_; + std::string path_; +}; + +void rethrow_thread_error(const std::exception_ptr& thread_error) +{ + if (thread_error) + { + std::rethrow_exception(thread_error); + } +} +} // namespace + +TEST(IpiSocketTest, WriteHeaderPadsToTwelveBytes) +{ + UnixSocketServer server; + std::string received; + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + char buffer[IPI_HEADER_LEN]; + recv_all(fd, buffer, sizeof(buffer)); + received.assign(buffer, sizeof(buffer)); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + socket.write_header("READY"); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); + EXPECT_EQ(padded_header("READY"), received); +} + +TEST(IpiSocketTest, CleanPeerCloseBeforeNextHeaderThrowsDedicatedSignal) +{ + UnixSocketServer server; + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + const std::string header = padded_header("STATUS"); + send_all(fd, header.data(), header.size()); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + EXPECT_EQ("STATUS", socket.read_header()); + EXPECT_THROW(socket.read_header(), IpiSocketClosed); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); +} + +TEST(IpiSocketTest, PartialHeaderCloseStaysRuntimeError) +{ + UnixSocketServer server; + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + const std::string partial = "STAT"; + send_all(fd, partial.data(), partial.size()); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + try + { + static_cast(socket.read_header()); + FAIL() << "partial header EOF should throw"; + } + catch (const IpiSocketClosed&) + { + FAIL() << "partial header EOF must not be treated as clean peer close"; + } + catch (const std::runtime_error& exc) + { + EXPECT_NE(std::string::npos, std::string(exc.what()).find("closed while reading header")); + } + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); +} From fe2c2a3cc61995b3383aaf7f2a6abf5ab1882907 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Thu, 9 Jul 2026 11:12:06 +0800 Subject: [PATCH 03/24] Fix socket driver CI regressions --- source/source_main/driver_run.cpp | 15 ++-- source/source_relax/socket_driver.cpp | 99 +++++++++++++++++++++++++-- 2 files changed, 101 insertions(+), 13 deletions(-) diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 3e7d669b6c2..408b507b902 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -38,9 +38,6 @@ void Driver::driver_run() { ModuleBase::TITLE("Driver", "driver_run"); - const std::string cal = PARAM.inp.calculation; - const bool socket_mode = (cal == "socket"); - //! 1: setup cell and atom information // this warning should not be here, mohan 2024-05-22 #ifndef __LCAO @@ -62,13 +59,17 @@ void Driver::driver_run() unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) - Input_para esolver_inp = PARAM.inp; - Input_para driver_inp = PARAM.inp; + const std::string cal = PARAM.inp.calculation; + const bool socket_mode = (cal == "socket"); + Input_para socket_esolver_inp = PARAM.inp; + Input_para socket_driver_inp = PARAM.inp; if (socket_mode) { - esolver_inp.calculation = "scf"; - driver_inp.calculation = cal; + socket_esolver_inp.calculation = "scf"; + socket_driver_inp.calculation = cal; } + const Input_para& esolver_inp = socket_mode ? socket_esolver_inp : PARAM.inp; + const Input_para& driver_inp = socket_mode ? socket_driver_inp : PARAM.inp; this->init_hardware(); diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index 02928fec9d3..bb0a6757f50 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -27,29 +27,52 @@ bool is_root() void bcast_double_vector(std::vector& values) { +#ifdef __MPI if (!values.empty()) { Parallel_Common::bcast_double(values.data(), static_cast(values.size())); } +#else + (void)values; +#endif +} + +void bcast_socket_int(int& value) +{ +#ifdef __MPI + Parallel_Common::bcast_int(value); +#else + (void)value; +#endif +} + +void bcast_socket_chars(char* value, const int size) +{ +#ifdef __MPI + Parallel_Common::bcast_char(value, size); +#else + (void)value; + (void)size; +#endif } void bcast_socket_string(std::string& value) { int size = static_cast(value.size()); - Parallel_Common::bcast_int(size); + bcast_socket_int(size); if (!is_root()) { value.resize(static_cast(size)); } if (size > 0) { - Parallel_Common::bcast_char(&value[0], size); + bcast_socket_chars(&value[0], size); } } void quit_if_root_io_failed(int root_failed, std::string root_message) { - Parallel_Common::bcast_int(root_failed); + bcast_socket_int(root_failed); bcast_socket_string(root_message); if (root_failed != 0) { @@ -85,6 +108,58 @@ std::vector ipi_cell_bohr_from_unitcell(const UnitCell& ucell) }; } +double max_wrapped_direct_delta_from_unitcell(const UnitCell& ucell, const std::vector& positions_bohr) +{ + if (positions_bohr.size() != static_cast(3 * ucell.nat)) + { + return 1.0e99; + } + + double out = 0.0; + int iat = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + const Atom* atom = &ucell.atoms[it]; + for (int ia = 0; ia < atom->na; ++ia) + { + const double tau_x = positions_bohr[3 * iat + 0] / ucell.lat0; + const double tau_y = positions_bohr[3 * iat + 1] / ucell.lat0; + const double tau_z = positions_bohr[3 * iat + 2] / ucell.lat0; + + double dx = 0.0; + double dy = 0.0; + double dz = 0.0; + ModuleBase::Mathzone::Cartesian_to_Direct(tau_x, + tau_y, + tau_z, + ucell.latvec.e11, + ucell.latvec.e12, + ucell.latvec.e13, + ucell.latvec.e21, + ucell.latvec.e22, + ucell.latvec.e23, + ucell.latvec.e31, + ucell.latvec.e32, + ucell.latvec.e33, + dx, + dy, + dz); + + double ddx = dx - atom->taud[ia].x; + double ddy = dy - atom->taud[ia].y; + double ddz = dz - atom->taud[ia].z; + ddx -= std::round(ddx); + ddy -= std::round(ddy); + ddz -= std::round(ddz); + out = std::max(out, std::abs(ddx)); + out = std::max(out, std::abs(ddy)); + out = std::max(out, std::abs(ddz)); + ++iat; + } + } + return out; +} + double max_abs_delta(const std::vector& a, const std::vector& b) { if (a.size() != b.size()) @@ -214,6 +289,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, std::vector virial_hartree(9, 0.0); const std::vector reference_cell = ipi_cell_bohr_from_unitcell(ucell); + bool checked_initial_positions = false; while (true) { @@ -306,8 +382,8 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } } quit_if_root_io_failed(io_failed, io_message); - Parallel_Common::bcast_int(rid); - Parallel_Common::bcast_int(nbytes); + bcast_socket_int(rid); + bcast_socket_int(nbytes); if (nbytes > 0 && is_root()) { ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; @@ -352,7 +428,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, quit_if_root_io_failed(io_failed, io_message); bcast_double_vector(cell); bcast_double_vector(inv_cell); - Parallel_Common::bcast_int(nat_socket); + bcast_socket_int(nat_socket); if (!is_root()) { positions.assign(static_cast(3 * nat_socket), 0.0); @@ -368,6 +444,17 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { ModuleBase::WARNING_QUIT("ABACUS socket", "variable-cell socket updates are not supported yet."); } + if (!checked_initial_positions) + { + checked_initial_positions = true; + if (max_wrapped_direct_delta_from_unitcell(ucell, positions) > 1.0e-5 && is_root()) + { + ModuleBase::WARNING( + "ABACUS socket", + "first POSDATA positions are not PBC-equivalent to STRU atom order; " + "i-PI POSDATA carries no species, so the client atoms should use the same atom order as STRU."); + } + } set_positions_from_ipi_bohr(ucell, positions); p_esolver->runner(ucell, istep); From 36353412e892d8d3a955d58f057285e886d1659d Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Thu, 9 Jul 2026 11:23:16 +0800 Subject: [PATCH 04/24] Reduce socket driver global dependencies --- source/source_main/driver_run.cpp | 9 +++------ source/source_relax/socket_driver.cpp | 12 ++++++++++-- source/source_relax/socket_driver.h | 12 ++++++++---- 3 files changed, 21 insertions(+), 12 deletions(-) diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index 408b507b902..e565871d604 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -59,17 +59,13 @@ void Driver::driver_run() unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) - const std::string cal = PARAM.inp.calculation; - const bool socket_mode = (cal == "socket"); Input_para socket_esolver_inp = PARAM.inp; - Input_para socket_driver_inp = PARAM.inp; + const bool socket_mode = (socket_esolver_inp.calculation == "socket"); if (socket_mode) { socket_esolver_inp.calculation = "scf"; - socket_driver_inp.calculation = cal; } const Input_para& esolver_inp = socket_mode ? socket_esolver_inp : PARAM.inp; - const Input_para& driver_inp = socket_mode ? socket_driver_inp : PARAM.inp; this->init_hardware(); @@ -84,6 +80,7 @@ void Driver::driver_run() #endif //! 4: different types of calculations + const std::string cal = PARAM.inp.calculation; if (cal == "md") { Run_MD::md_line(ucell, p_esolver, PARAM); @@ -91,7 +88,7 @@ void Driver::driver_run() else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf" || cal == "socket") { Relax_Driver rl_driver; - rl_driver.relax_driver(p_esolver, ucell, driver_inp, GlobalV::ofs_running); + rl_driver.relax_driver(p_esolver, ucell, PARAM.inp, GlobalV::ofs_running); } else if (cal == "get_s") { diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index bb0a6757f50..ecd4d1be27e 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -2,11 +2,13 @@ #include "source_relax/ipi_socket.h" #include "source_base/global_function.h" -#include "source_base/global_variable.h" #include "source_base/mathzone.h" #include "source_base/parallel_common.h" #include "source_base/timer.h" +#include "source_cell/unitcell.h" #include "source_cell/update_cell.h" +#include "source_esolver/esolver.h" +#include "source_io/module_parameter/input_parameter.h" #include #include @@ -22,7 +24,13 @@ constexpr int IPI_RANK_ROOT = 0; bool is_root() { - return GlobalV::MY_RANK == IPI_RANK_ROOT; +#ifdef __MPI + int rank = IPI_RANK_ROOT; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + return rank == IPI_RANK_ROOT; +#else + return true; +#endif } void bcast_double_vector(std::vector& values) diff --git a/source/source_relax/socket_driver.h b/source/source_relax/socket_driver.h index 51b75289cb3..86c180fc42e 100644 --- a/source/source_relax/socket_driver.h +++ b/source/source_relax/socket_driver.h @@ -1,11 +1,15 @@ #ifndef ABACUS_SOURCE_RELAX_SOCKET_DRIVER_H #define ABACUS_SOURCE_RELAX_SOCKET_DRIVER_H -#include "source_cell/unitcell.h" -#include "source_esolver/esolver.h" -#include "source_io/module_parameter/input_parameter.h" +#include -#include +class UnitCell; +struct Input_para; + +namespace ModuleESolver +{ +class ESolver; +} class Socket_Driver { From 5711ab3ceccf6b7df036ca5cf6a932de38eb22a3 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 19 Jul 2026 13:20:52 +0800 Subject: [PATCH 05/24] Address socket i-PI review feedback --- docs/advanced/input_files/input-main.md | 13 +- docs/advanced/interface/ase.md | 58 +++++ docs/parameters.yaml | 13 +- interfaces/ASE_interface/README.md | 1 + interfaces/ASE_interface/abacuslite/core.py | 216 ++++++++++++++++++ interfaces/ASE_interface/examples/socketio.py | 157 +++++++++++++ .../module_parameter/input_parameter.h | 1 + .../read_input_item_system.cpp | 30 ++- .../test_serial/read_input_item_test.cpp | 33 ++- source/source_main/driver_run.cpp | 14 +- source/source_relax/CMakeLists.txt | 2 +- source/source_relax/relax_driver.cpp | 2 +- source/source_relax/socket_driver.cpp | 10 +- .../{ipi_socket.cpp => socket_ipi.cpp} | 6 +- .../{ipi_socket.h => socket_ipi.h} | 4 +- source/source_relax/test/CMakeLists.txt | 4 +- ...pi_socket_test.cpp => socket_ipi_test.cpp} | 4 +- 17 files changed, 531 insertions(+), 37 deletions(-) create mode 100644 interfaces/ASE_interface/examples/socketio.py rename source/source_relax/{ipi_socket.cpp => socket_ipi.cpp} (97%) rename source/source_relax/{ipi_socket.h => socket_ipi.h} (94%) rename source/source_relax/test/{ipi_socket_test.cpp => socket_ipi_test.cpp} (98%) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 529e2657edb..ee1fb6c5b32 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -9,6 +9,7 @@ - [suffix](#suffix) - [ntype](#ntype) - [calculation](#calculation) + - [socket\_driver](#socket_driver) - [esolver\_type](#esolver_type) - [symmetry](#symmetry) - [symmetry\_prec](#symmetry_prec) @@ -583,7 +584,6 @@ - relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations - cell-relax: perform cell relaxation calculations - md: perform molecular dynamics simulations - - socket: run as a socket client for external drivers using the i-PI protocol - get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information - get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information - get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -593,6 +593,15 @@ - test_neighbour: obtain information of neighboring atoms (for LCAO basis only), please specify a positive search_radius manually - **Default**: scf +### socket_driver + +- **Type**: Boolean +- **Availability**: *calculation==scf* +- **Description**: If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. + + > Note: Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used. +- **Default**: False + ### esolver_type - **Type**: String @@ -844,7 +853,7 @@ - **Description**: Charge extrapolation method for MD, relaxation, and socket-driven calculations. When set to default, ABACUS chooses second-order for md, first-order for - relax/cell-relax/socket, and atomic for other calculations. Socket-driven + relax/cell-relax and socket_driver calculations, and atomic for other calculations. Socket-driven molecular dynamics can explicitly set second-order if the external driver updates structures smoothly enough for second-order extrapolation. - **Default**: default diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index e9b7b062809..82289cbac33 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -103,6 +103,64 @@ In the new implementation, we limit the range of functionalties supported to mai Please read the examples in `interfaces/ASE_interface/examples/` for more details. +### Socket I/O with ASE + +For socket-driven ASE workflows, use the `AbacusSocketIO` calculator. ASE runs the i-PI socket server, while ABACUS keeps `calculation=scf` and is launched with `socket_driver=1` as the client. The protocol is simple: ASE sends atomic positions and cell data to ABACUS; ABACUS evaluates one SCF step for that structure and returns energy, forces, and virial. See the [ASE socket I/O documentation](https://ase-lib.org/ase/calculators/socketio/socketio.html) and the i-PI reference paper, [Ceriotti et al., Comput. Phys. Commun. 185, 1019-1026 (2014)](https://doi.org/10.1016/j.cpc.2013.10.027), for the protocol background. + +Build ABACUS as usual before using this interface. PW-only builds work with `basis_type=pw`; LCAO socket calculations require an LCAO-enabled executable. No extra socket library is required. + +With CMake, choose the executable according to the basis: + +```bash +cmake -S . -B build-pw -DENABLE_MPI=ON -DENABLE_LCAO=OFF +cmake --build build-pw --target abacus_pw_para -j + +cmake -S . -B build-lcao -DENABLE_MPI=ON -DENABLE_LCAO=ON +cmake --build build-lcao --target abacus_basic_para -j +``` + +With the ABACUS toolchain workflow, build the normal ABACUS executable with LCAO support when `basis_type=lcao` is needed, then pass that executable to `AbacusProfile(command=...)`. The ASE interface can be installed from this repository with: + +```bash +cd interfaces/ASE_interface +pip install . +``` + +A minimal socket calculator setup is: + +```python +from ase.optimize import BFGS +from abacuslite import AbacusProfile, AbacusSocketIO + +aprof = AbacusProfile( + command="mpirun -np 4 /path/to/abacus", + pseudo_dir="/path/to/pseudopotentials", + orbital_dir="/path/to/orbitals", + omp_num_threads=1, +) + +abacus = AbacusSocketIO( + profile=aprof, + directory="socketio", + unixsocket="abacus_si", + pseudopotentials={"Si": "Si_ONCV_PBE-1.0.upf"}, + basissets={"Si": "Si_gga_8au_100Ry_2s2p1d.orb"}, + inp={"calculation": "scf", "basis_type": "lcao", "kspacing": 0.1}, +) + +with abacus as calc: + atoms.calc = calc + BFGS(atoms).run(fmax=0.05) +``` + +`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. The socket endpoint is selected by the calculator arguments and passed to ABACUS through `ABACUS_SOCKET_ADDRESS`, for example `localhost:31415` or `/tmp/ipi_abacus_si:UNIX`. Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. + +A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. + +The i-PI protocol does not transmit element symbols. `AbacusSocketIO` therefore sorts the internal socket atoms with the same first-occurrence species grouping used when writing `STRU`, and maps returned forces back to the original ASE `Atoms` order. This avoids silent force/atom mismatches when structures are read from CIF, extxyz, POSCAR, or other formats whose atom order is not already grouped for ABACUS. Users should not manually reorder atoms for socket I/O; pass the physical ASE `Atoms` object directly to the calculator. + +A complete fixed-cell validation and benchmark example is available in `interfaces/ASE_interface/examples/socketio.py`. + ## SPAP Analysis [SPAP](https://github.com/chuanxun/StructurePrototypeAnalysisPackage) (Structure Prototype Analysis Package) is written by Dr. Chuanxun Su to analyze symmetry and compare similarity of large amount of atomic structures. The coordination characterization function (CCF) is used to diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 2298e10d111..9c9a250707b 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -28,7 +28,6 @@ parameters: * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations * cell-relax: perform cell relaxation calculations * md: perform molecular dynamics simulations - * socket: run as a socket client for external drivers using the i-PI protocol * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -39,6 +38,16 @@ parameters: default_value: scf unit: "" availability: "" + - name: socket_driver + category: System variables + type: Boolean + description: | + If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. + + [NOTE] Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used. + default_value: "False" + unit: "" + availability: calculation==scf - name: esolver_type category: System variables type: String @@ -333,7 +342,7 @@ parameters: Charge extrapolation method for MD, relaxation, and socket-driven calculations. When set to default, ABACUS chooses second-order for md, first-order for - relax/cell-relax/socket, and atomic for other calculations. Socket-driven + relax/cell-relax and socket_driver calculations, and atomic for other calculations. Socket-driven molecular dynamics can explicitly set second-order if the external driver updates structures smoothly enough for second-order extrapolation. default_value: default diff --git a/interfaces/ASE_interface/README.md b/interfaces/ASE_interface/README.md index 824e69900de..a76e68e177e 100644 --- a/interfaces/ASE_interface/README.md +++ b/interfaces/ASE_interface/README.md @@ -32,6 +32,7 @@ Please refer to the example scripts in the `examples` folder. Recommended learni 7. **constraintmd.py** - Constrained molecular dynamics simulation 8. **metadynamics.py** - Metadynamics simulation 9. **neb.py** - Nudged Elastic Band (NEB) calculation +10. **socketio.py** - ASE optimization with the AbacusSocketIO calculator running ABACUS as an i-PI socket client More usage examples will be provided in future versions. diff --git a/interfaces/ASE_interface/abacuslite/core.py b/interfaces/ASE_interface/abacuslite/core.py index a1aad1e773b..e2c00d4c889 100644 --- a/interfaces/ASE_interface/abacuslite/core.py +++ b/interfaces/ASE_interface/abacuslite/core.py @@ -45,6 +45,7 @@ GenericFileIOCalculator, read_stdout ) +from ase.calculators.socketio import SocketIOCalculator from ase.atoms import Atoms from ase.dft.kpoints import BandPath from ase.io import read @@ -135,6 +136,21 @@ def get_calculator_command(self, inputfile) -> List[str]: # additional inputfile argument is not used. return [] + def socketio_argv_inet(self, port: Optional[int] = None) -> List[str]: + port = 31415 if port is None else port + return [ + 'env', + f'ABACUS_SOCKET_ADDRESS=localhost:{port}', + *self._split_command, + ] + + def socketio_argv_unix(self, socket: str) -> List[str]: + return [ + 'env', + f'ABACUS_SOCKET_ADDRESS=/tmp/ipi_{socket}:UNIX', + *self._split_command, + ] + def version(self) -> str: '''get the abacus version information''' cmd_ = [*self._split_command, '--version'] @@ -443,6 +459,17 @@ def __init__(self, directory=directory, ) + def write_input(self, atoms, properties=None, system_changes=None): + if properties is None: + properties = self.template.implemented_properties + self.template.write_input( + profile=self.profile, + directory=Path(self.directory), + atoms=atoms, + parameters=self.parameters, + properties=properties, + ) + @classmethod def restart(cls, profile=None, directory='.', **kwargs): '''instantiate one ABACUS calculator from an existing job directory, @@ -557,11 +584,200 @@ def band_structure(self, efermi=None): from ase.spectrum.band_structure import get_band_structure return get_band_structure(calc=self, reference=efermi) +class AbacusSocketIO(SocketIOCalculator): + """ASE socket I/O calculator that launches ABACUS as an i-PI client. + + A socket calculator owns one ABACUS process with one fixed INPUT/STRU + setup. The i-PI protocol can update positions, but electronic-structure + parameters such as k-points, spin, basis, pseudopotentials, and species + require a new calculator instance. Energy-only ASE calls are accepted, but + ABACUS still computes forces because i-PI GETFORCE returns energy, forces, + and virial together. + """ + + def __init__(self, + profile=None, + directory='.', + port=None, + unixsocket=None, + timeout=None, + log=None, + **kwargs): + inp = self._socket_inp(kwargs.pop('inp', {})) + self.abacus = Abacus( + profile=profile, + directory=directory, + inp=inp, + **kwargs, + ) + self._reference_cell = None + super().__init__( + port=port, + unixsocket=unixsocket, + timeout=timeout, + log=log, + launch_client=self._launch_client, + ) + + def calculate(self, atoms=None, properties=['energy'], system_changes=None): + from ase.calculators.calculator import ( + PropertyNotImplementedError, + all_changes, + ) + from ase.stress import full_3x3_to_voigt_6_stress + + if system_changes is None: + system_changes = all_changes + if atoms is None: + atoms = self.atoms + if atoms is None: + raise ValueError('AbacusSocketIO.calculate requires atoms') + + bad = [change for change in system_changes + if change not in self.supported_changes] + if self.atoms is not None and any(bad): + raise PropertyNotImplementedError( + 'Cannot change {} through IPI protocol. ' + 'Please create new socket calculator.' + .format(bad if len(bad) > 1 else bad[0])) + + self._check_fixed_cell(atoms) + order = self._socket_sort_indices(atoms) + socket_atoms = atoms[order] + self.atoms = atoms.copy() + + if self.server is None: + self.server = self.launch_server() + proc = self.launch_client(socket_atoms, properties, + port=self._port, + unixsocket=self._unixsocket) + self.server.proc = proc + + results = self.server.calculate(socket_atoms) + results['free_energy'] = results['energy'] + virial = results.pop('virial') + if self.atoms.cell.rank == 3 and any(self.atoms.pbc): + vol = atoms.get_volume() + results['stress'] = -full_3x3_to_voigt_6_stress(virial) / vol + if 'forces' in results: + results['forces'] = self._forces_to_input_order( + results['forces'], order) + self.results.update(results) + + def _check_fixed_cell(self, atoms): + from ase.calculators.calculator import PropertyNotImplementedError + + cell = atoms.cell.array.copy() + if self._reference_cell is None: + self._reference_cell = cell + return + max_delta = np.max(np.abs(cell - self._reference_cell)) + if max_delta > 1.0e-10: + raise PropertyNotImplementedError( + 'AbacusSocketIO is fixed-cell only; create a new socket ' + 'calculator for a changed cell, or use the normal Abacus ' + 'FileIO calculator for variable-cell workflows.' + ) + + def set(self, **kwargs): + if kwargs: + raise ValueError( + 'AbacusSocketIO input parameters are fixed after construction; ' + 'create a new AbacusSocketIO calculator to change k-points, ' + 'spin, basis, pseudopotentials, species, or other INPUT/STRU ' + 'settings.' + ) + return super().set(**kwargs) + + def _launch_client(self, atoms, properties=None, port=None, unixsocket=None): + from subprocess import Popen + + if properties is None: + properties = self.abacus.template.implemented_properties + + directory = Path(self.abacus.directory) + directory.mkdir(exist_ok=True, parents=True) + + if hasattr(self.abacus, 'write_inputfiles'): + self.abacus.write_inputfiles(atoms, properties) + else: + self.abacus.write_input(atoms, properties=properties) + + if unixsocket is not None: + argv = self.abacus.profile.socketio_argv_unix(socket=unixsocket) + else: + argv = self.abacus.profile.socketio_argv_inet(port=port) + + stdout = open(directory / self.abacus.template.outputname, 'w') + stderr = open(directory / self.abacus.template.errorname, 'w') + try: + return Popen(argv, cwd=directory, env=os.environ, + stdout=stdout, stderr=stderr) + finally: + stdout.close() + stderr.close() + + @staticmethod + def _socket_inp(inp): + inp = dict(inp) + calculation = inp.get('calculation', 'scf') + if calculation != 'scf': + raise ValueError('ABACUS socket I/O requires calculation="scf"') + inp.update({ + 'calculation': 'scf', + 'socket_driver': 1, + 'cal_force': 1, + }) + return inp + + @staticmethod + def _socket_sort_indices(atoms): + return species_group_indices(atoms.get_chemical_symbols()) + + @staticmethod + def _forces_to_input_order(forces, order): + reordered = np.empty_like(forces) + for sorted_index, original_index in enumerate(order): + reordered[original_index] = forces[sorted_index] + return reordered + + class TestAbacusCalculator(unittest.TestCase): here = Path(__file__).parent pporb = here.parent.parent.parent / 'tests' / 'PP_ORB' + def test_socketio_species_order_mapping(self): + atoms = Atoms(symbols=['Si', 'O', 'C', 'Si', 'O', 'C']) + order = AbacusSocketIO._socket_sort_indices(atoms) + self.assertEqual(order, [0, 3, 1, 4, 2, 5]) + + socket_forces = np.arange(18).reshape(6, 3) + input_forces = AbacusSocketIO._forces_to_input_order( + socket_forces, order) + + expected = np.empty_like(socket_forces) + for sorted_index, original_index in enumerate(order): + expected[original_index] = socket_forces[sorted_index] + np.testing.assert_array_equal(input_forces, expected) + + def test_socketio_rejects_parameter_changes(self): + calc = object.__new__(AbacusSocketIO) + with self.assertRaisesRegex(ValueError, 'fixed after construction'): + calc.set(kpts={'mode': 'mp-sampling', 'nk': [2, 2, 2]}) + + def test_socketio_rejects_cell_changes(self): + from ase.calculators.calculator import PropertyNotImplementedError + + calc = object.__new__(AbacusSocketIO) + calc.atoms = Atoms('Si', cell=[5.0, 5.0, 5.0], pbc=True) + calc._reference_cell = calc.atoms.cell.array.copy() + + changed = calc.atoms.copy() + changed.cell[0, 0] = 5.1 + with self.assertRaisesRegex(PropertyNotImplementedError, 'fixed-cell'): + calc._check_fixed_cell(changed) + def test_calculator_results(self): from ase.build.bulk import bulk silicon = bulk('Si', crystalstructure='diamond', a=5.43) diff --git a/interfaces/ASE_interface/examples/socketio.py b/interfaces/ASE_interface/examples/socketio.py new file mode 100644 index 00000000000..e3087453ffc --- /dev/null +++ b/interfaces/ASE_interface/examples/socketio.py @@ -0,0 +1,157 @@ +""" +This example validates and benchmarks ABACUS socket I/O from ASE. + +ASE runs as the i-PI socket server and ABACUS runs as the socket client. +ABACUS keeps calculation=scf and enables socket_driver internally. + +The script checks two PR-review relevant points: +1. socket SCF gives the same energy and forces as a normal non-socket SCF; +2. repeated socket calculations avoid relaunching ABACUS and are faster than + the normal FileIO calculator for a sequence of SCF force evaluations. + +The i-PI protocol does not carry element symbols. AbacusSocketIO handles +the required STRU/socket atom-order alignment internally and returns forces +in the original ASE Atoms order. +""" +import os +import shutil +import time +from pathlib import Path + +import numpy as np +from ase import Atoms +from abacuslite import Abacus, AbacusProfile, AbacusSocketIO + +here = Path(__file__).parent +pporb = here.parent.parent.parent / 'tests' / 'PP_ORB' + +aprof = AbacusProfile( + command=os.environ.get('ABACUS_COMMAND', 'mpirun -np 4 abacus'), + pseudo_dir=pporb, + orbital_dir=pporb, + omp_num_threads=1, +) + +common_kwargs = { + 'pseudopotentials': {'Si': 'Si_ONCV_PBE-1.0.upf'}, + 'basissets': {'Si': 'Si_gga_8au_100Ry_2s2p1d.orb'}, + 'inp': { + 'calculation': 'scf', + 'nspin': 1, + 'basis_type': 'lcao', + 'ks_solver': 'scalapack_gvx', + 'ecutwfc': 30, + 'symmetry': 0, + 'kspacing': 0.5, + 'scf_thr': 1e-8, + 'scf_nmax': 40, + 'chg_extrap': 'atomic', + 'cal_force': 1, + }, +} + +base_atoms = Atoms( + 'Si2', + positions=[[0.0, 0.0, 0.0], [1.25, 1.25, 1.25]], + cell=[5.43, 5.43, 5.43], + pbc=True, +) + + +def clean(directory): + shutil.rmtree(directory, ignore_errors=True) + + +def run_fileio(atoms, directory): + clean(directory) + calc = Abacus(profile=aprof, directory=str(directory), **common_kwargs) + atoms = atoms.copy() + atoms.calc = calc + forces = atoms.get_forces() + energy = atoms.get_potential_energy() + return energy, forces + + +def run_socketio(atoms, directory, socket_name): + clean(directory) + calc = AbacusSocketIO( + profile=aprof, + directory=str(directory), + unixsocket=socket_name, + timeout=120, + **common_kwargs, + ) + atoms = atoms.copy() + with calc: + atoms.calc = calc + energy = atoms.get_potential_energy() + forces = atoms.get_forces() + return energy, forces + + +def displaced_structures(): + structures = [] + for scale in (0.00, 0.03, -0.02, 0.05): + atoms = base_atoms.copy() + atoms.positions[1] += scale + structures.append(atoms) + return structures + + +fileio_dir = here / 'socketio_fileio_scf' +socket_dir = here / 'socketio_socket_scf' +bench_fileio_dir = here / 'socketio_bench_fileio' +bench_socket_dir = here / 'socketio_bench_socket' + +try: + reference_energy, reference_forces = run_fileio(base_atoms, fileio_dir) + socket_energy, socket_forces = run_socketio( + base_atoms, socket_dir, 'abacus_si_check') + + energy_diff = abs(socket_energy - reference_energy) + force_diff = np.max(np.abs(socket_forces - reference_forces)) + print(f'FileIO SCF energy: {reference_energy:.12f} eV') + print(f'Socket SCF energy: {socket_energy:.12f} eV') + print(f'|dE|: {energy_diff:.3e} eV') + print(f'max |dF|: {force_diff:.3e} eV/Angstrom') + assert energy_diff < 1e-4 + assert force_diff < 1e-5 + + structures = displaced_structures() + + clean(bench_fileio_dir) + fileio_calc = Abacus( + profile=aprof, + directory=str(bench_fileio_dir), + **common_kwargs, + ) + t0 = time.perf_counter() + for atoms in structures: + atoms = atoms.copy() + atoms.calc = fileio_calc + atoms.get_forces() + fileio_seconds = time.perf_counter() - t0 + + clean(bench_socket_dir) + socket_calc = AbacusSocketIO( + profile=aprof, + directory=str(bench_socket_dir), + unixsocket='abacus_si_bench', + timeout=120, + **common_kwargs, + ) + t0 = time.perf_counter() + with socket_calc as calc: + for atoms in structures: + atoms = atoms.copy() + atoms.calc = calc + atoms.get_forces() + socket_seconds = time.perf_counter() - t0 + + speedup = fileio_seconds / socket_seconds + print(f'FileIO repeated SCF force time: {fileio_seconds:.2f} s') + print(f'Socket repeated SCF force time: {socket_seconds:.2f} s') + print(f'Socket speedup vs FileIO: {speedup:.2f}x') +finally: + for directory in (fileio_dir, socket_dir, bench_fileio_dir, bench_socket_dir): + clean(directory) diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 3045d663361..c7cb990da21 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -19,6 +19,7 @@ struct Input_para std::string calculation = "scf"; ///< "scf" : self consistent calculation. ///< "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 /* symmetry level: -1, no symmetry at all; diff --git a/source/source_io/module_parameter/read_input_item_system.cpp b/source/source_io/module_parameter/read_input_item_system.cpp index 6ed572b21f9..bc2d1c94d4f 100644 --- a/source/source_io/module_parameter/read_input_item_system.cpp +++ b/source/source_io/module_parameter/read_input_item_system.cpp @@ -75,7 +75,7 @@ void ReadInput::item_system() } { Input_Item item("calculation"); - item.annotation = "scf; relax; md; socket; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; + item.annotation = "scf; relax; md; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour"; item.category = "System variables"; item.type = "String"; item.description = R"(Specify the type of calculation. @@ -85,7 +85,6 @@ void ReadInput::item_system() * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations * cell-relax: perform cell relaxation calculations * md: perform molecular dynamics simulations -* socket: run as a socket client for external drivers using the i-PI protocol * get_pchg: obtain partial (band-decomposed) charge densities (for LCAO basis only). See out_pchg for more information * get_wf: obtain real space wave functions (for LCAO basis only). See out_wfc_norm and out_wfc_re_im for more information * get_s: obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is SR.csr with file format being the same as that generated by out_mat_hs2 @@ -103,7 +102,6 @@ void ReadInput::item_system() std::vector callist = {"scf", "relax", "md", - "socket", "cell-relax", "nscf", "get_s", @@ -138,6 +136,24 @@ void ReadInput::item_system() sync_string(input.calculation); this->add_item(item); } + { + Input_Item item("socket_driver"); + item.annotation = "run as a socket client for external drivers using the i-PI protocol"; + item.category = "System variables"; + item.type = "Boolean"; + item.description = R"(If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. + +[NOTE] Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used.)"; + item.default_value = "False"; + read_sync_bool(input.socket_driver); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.socket_driver && para.input.calculation != "scf") + { + ModuleBase::WARNING_QUIT("ReadInput", "socket_driver is only supported with calculation = scf."); + } + }; + this->add_item(item); + } { Input_Item item("esolver_type"); item.annotation = "the energy solver: ksdft, sdft, ofdft, tdofdft, tddft, lj, dp, ks-lr, lr"; @@ -268,9 +284,9 @@ void ReadInput::item_system() item.description = "If set to True, calculate the force at the end of the electronic iteration."; item.default_value = "False"; item.reset_value = [](const Input_Item& item, Parameter& para) { - std::vector use_force = {"cell-relax", "relax", "md", "socket"}; + std::vector use_force = {"cell-relax", "relax", "md"}; std::vector not_use_force = {"get_wf", "get_pchg", "get_s"}; - if (std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) + if (para.input.socket_driver || std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) { if (!para.input.cal_force) { @@ -878,7 +894,7 @@ Available options are: item.description = R"(Charge extrapolation method for MD, relaxation, and socket-driven calculations. When set to default, ABACUS chooses second-order for md, first-order for -relax/cell-relax/socket, and atomic for other calculations. Socket-driven +relax/cell-relax and socket_driver calculations, and atomic for other calculations. Socket-driven molecular dynamics can explicitly set second-order if the external driver updates structures smoothly enough for second-order extrapolation.)"; item.default_value = "default"; @@ -889,7 +905,7 @@ updates structures smoothly enough for second-order extrapolation.)"; para.input.chg_extrap = "second-order"; } else if (para.input.chg_extrap == "default" - && (para.input.calculation == "relax" || para.input.calculation == "cell-relax" || para.input.calculation == "socket")) + && (para.input.calculation == "relax" || para.input.calculation == "cell-relax" || para.input.socket_driver)) { para.input.chg_extrap = "first-order"; } diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index eb043fe1b03..0a53fa13d85 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -64,6 +64,27 @@ TEST_F(InputTest, Item_test) EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + + param.input.calculation = "socket"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + } + + { // socket_driver + auto it = find_label("socket_driver", readinput.input_lists); + param.input.socket_driver = true; + param.input.calculation = "nscf"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NOTICE")); + + param.input.socket_driver = true; + param.input.calculation = "scf"; + EXPECT_NO_THROW(it->second.check_value(it->second, param)); + param.input.socket_driver = false; } { // esolver_type @@ -246,6 +267,7 @@ TEST_F(InputTest, Item_test) auto it = find_label("cal_force", readinput.input_lists); param.input.calculation = "cell-relax"; param.input.cal_force = false; + param.input.socket_driver = false; it->second.reset_value(it->second, param); EXPECT_EQ(param.input.cal_force, true); @@ -253,6 +275,13 @@ TEST_F(InputTest, Item_test) param.input.cal_force = true; it->second.reset_value(it->second, param); EXPECT_EQ(param.input.cal_force, false); + + param.input.calculation = "scf"; + param.input.socket_driver = true; + param.input.cal_force = false; + it->second.reset_value(it->second, param); + EXPECT_EQ(param.input.cal_force, true); + param.input.socket_driver = false; } { // ecutrho auto it = find_label("ecutrho", readinput.input_lists); @@ -363,12 +392,14 @@ TEST_F(InputTest, Item_test) EXPECT_EQ(param.input.chg_extrap, "first-order"); param.input.chg_extrap = "default"; - param.input.calculation = "socket"; + param.input.calculation = "scf"; + param.input.socket_driver = true; it->second.reset_value(it->second, param); EXPECT_EQ(param.input.chg_extrap, "first-order"); param.input.chg_extrap = "default"; param.input.calculation = "none"; + param.input.socket_driver = false; it->second.reset_value(it->second, param); EXPECT_EQ(param.input.chg_extrap, "atomic"); diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index e565871d604..d78ea0058e8 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -59,20 +59,12 @@ void Driver::driver_run() unitcell::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) - Input_para socket_esolver_inp = PARAM.inp; - const bool socket_mode = (socket_esolver_inp.calculation == "socket"); - if (socket_mode) - { - socket_esolver_inp.calculation = "scf"; - } - const Input_para& esolver_inp = socket_mode ? socket_esolver_inp : PARAM.inp; - this->init_hardware(); - ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(esolver_inp, ucell); + ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); //! 3: initialize Esolver and fill json-structure - p_esolver->before_all_runners(ucell, esolver_inp); + p_esolver->before_all_runners(ucell, PARAM.inp); // this Json part should be moved to before_all_runners, mohan 2024-05-12 #ifdef __RAPIDJSON @@ -85,7 +77,7 @@ void Driver::driver_run() { Run_MD::md_line(ucell, p_esolver, PARAM); } - else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf" || cal == "socket") + else if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf") { Relax_Driver rl_driver; rl_driver.relax_driver(p_esolver, ucell, PARAM.inp, GlobalV::ofs_running); diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 84f69c06c80..77079eb6a1f 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -2,7 +2,7 @@ add_library( relax OBJECT relax_data.cpp - ipi_socket.cpp + socket_ipi.cpp socket_driver.cpp cg_base.cpp relax_driver.cpp diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 0cada1bfc3e..922a230d882 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -18,7 +18,7 @@ void Relax_Driver::relax_driver( ModuleBase::TITLE("Relax_Driver", "relax_driver"); ModuleBase::timer::start("Relax_Driver", "relax_driver"); - if (inp.calculation == "socket") + if (inp.socket_driver) { Socket_Driver socket_driver; socket_driver.socket_driver(p_esolver, ucell, inp, ofs_running); diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index ecd4d1be27e..42fafda5e63 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -1,6 +1,6 @@ #include "socket_driver.h" -#include "source_relax/ipi_socket.h" +#include "source_relax/socket_ipi.h" #include "source_base/global_function.h" #include "source_base/mathzone.h" #include "source_base/parallel_common.h" @@ -94,9 +94,9 @@ std::string bcast_header(std::string header) return header; } -std::string ipi_address() +std::string socket_address() { - const char* env = std::getenv("ABACUS_IPI_ADDRESS"); + const char* env = std::getenv("ABACUS_SOCKET_ADDRESS"); if (env == nullptr || std::string(env).empty()) { return "localhost:31415"; @@ -263,7 +263,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } if (!inp.cal_force) { - ModuleBase::WARNING_QUIT("ABACUS socket", "socket calculation requires cal_force=1 for i-PI GETFORCE."); + ModuleBase::WARNING_QUIT("ABACUS socket", "socket_driver requires cal_force=1 for i-PI GETFORCE."); } IpiSocket socket; @@ -276,7 +276,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { try { - const std::string address = ipi_address(); + const std::string address = socket_address(); ofs_running << " ABACUS socket driver connecting to i-PI endpoint " << address << std::endl; socket.connect(address); } diff --git a/source/source_relax/ipi_socket.cpp b/source/source_relax/socket_ipi.cpp similarity index 97% rename from source/source_relax/ipi_socket.cpp rename to source/source_relax/socket_ipi.cpp index 304a0181813..2bb09385b0d 100644 --- a/source/source_relax/ipi_socket.cpp +++ b/source/source_relax/socket_ipi.cpp @@ -1,4 +1,4 @@ -#include "source_relax/ipi_socket.h" +#include "source_relax/socket_ipi.h" #include #include @@ -149,6 +149,10 @@ std::string IpiSocket::read_header() { continue; } + if (errno == ECONNRESET && done == 0) + { + throw IpiSocketClosed("i-PI socket peer reset before next header"); + } throw std::runtime_error(errno_message("i-PI socket header read failed")); } done += static_cast(nread); diff --git a/source/source_relax/ipi_socket.h b/source/source_relax/socket_ipi.h similarity index 94% rename from source/source_relax/ipi_socket.h rename to source/source_relax/socket_ipi.h index 28ce2cb81ff..eac15eb41ac 100644 --- a/source/source_relax/ipi_socket.h +++ b/source/source_relax/socket_ipi.h @@ -1,5 +1,5 @@ -#ifndef ABACUS_IPI_SOCKET_H -#define ABACUS_IPI_SOCKET_H +#ifndef ABACUS_SOCKET_IPI_H +#define ABACUS_SOCKET_IPI_H #include #include diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index d1c867dc690..51dc9a82d3b 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -8,8 +8,8 @@ install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( - TARGET MODULE_RELAX_ipi_socket_test - SOURCES ipi_socket_test.cpp ../ipi_socket.cpp + TARGET MODULE_RELAX_socket_ipi_test + SOURCES socket_ipi_test.cpp ../socket_ipi.cpp ) AddTest( diff --git a/source/source_relax/test/ipi_socket_test.cpp b/source/source_relax/test/socket_ipi_test.cpp similarity index 98% rename from source/source_relax/test/ipi_socket_test.cpp rename to source/source_relax/test/socket_ipi_test.cpp index 21db9be2b90..930ab1fe629 100644 --- a/source/source_relax/test/ipi_socket_test.cpp +++ b/source/source_relax/test/socket_ipi_test.cpp @@ -1,4 +1,4 @@ -#include "../ipi_socket.h" +#include "../socket_ipi.h" #include "gtest/gtest.h" @@ -80,7 +80,7 @@ class UnixSocketServer public: UnixSocketServer() { - char dir_template[] = "/tmp/abacus_ipi_socket_test_XXXXXX"; + char dir_template[] = "/tmp/abacus_socket_ipi_test_XXXXXX"; char* made_dir = ::mkdtemp(dir_template); if (made_dir == nullptr) { From faafda3855a1bd7af3ebba43c33fc79bb078d961 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 19 Jul 2026 14:00:19 +0800 Subject: [PATCH 06/24] Fix socket parameter docs sync --- docs/advanced/interface/ase.md | 2 +- docs/parameters.yaml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index 82289cbac33..9e4153e3572 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -153,7 +153,7 @@ with abacus as calc: BFGS(atoms).run(fmax=0.05) ``` -`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. The socket endpoint is selected by the calculator arguments and passed to ABACUS through `ABACUS_SOCKET_ADDRESS`, for example `localhost:31415` or `/tmp/ipi_abacus_si:UNIX`. Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. +`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. The socket endpoint is selected by the calculator arguments and passed to ABACUS through `ABACUS_SOCKET_ADDRESS`. `unixsocket="abacus_si"` means that ASE listens on the local Unix-domain socket `/tmp/ipi_abacus_si`; abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=/tmp/ipi_abacus_si:UNIX`. For TCP sockets, use the `port` argument instead, which maps to an address such as `localhost:31415`. Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 9c9a250707b..74127d1b841 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -47,7 +47,7 @@ parameters: [NOTE] Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used. default_value: "False" unit: "" - availability: calculation==scf + availability: "" - name: esolver_type category: System variables type: String From 9f26e9f40188afc082c7e0937f702e416dea7c37 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 19 Jul 2026 14:06:28 +0800 Subject: [PATCH 07/24] Clarify socket endpoint address formats --- docs/advanced/input_files/input-main.md | 7 ++++++- docs/advanced/interface/ase.md | 11 +++++++++-- docs/parameters.yaml | 5 ++++- .../module_parameter/read_input_item_system.cpp | 5 ++++- 4 files changed, 23 insertions(+), 5 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index ee1fb6c5b32..dee59db41d1 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -599,7 +599,12 @@ - **Availability**: *calculation==scf* - **Description**: If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. - > Note: Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used. + > Note: Use calculation = scf with socket_driver = True. ABACUS connects to the external i-PI server selected by ABACUS_SOCKET_ADDRESS. If ABACUS_SOCKET_ADDRESS is unset, ABACUS uses localhost:31415. The value can use one of two forms: + + - host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. + - path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. + + When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. - **Default**: False ### esolver_type diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index 9e4153e3572..7e4ed1a6af5 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -119,7 +119,7 @@ cmake -S . -B build-lcao -DENABLE_MPI=ON -DENABLE_LCAO=ON cmake --build build-lcao --target abacus_basic_para -j ``` -With the ABACUS toolchain workflow, build the normal ABACUS executable with LCAO support when `basis_type=lcao` is needed, then pass that executable to `AbacusProfile(command=...)`. The ASE interface can be installed from this repository with: +With the ABACUS toolchain workflow, build the normal ABACUS executable with LCAO support when `basis_type=lcao` is needed, then pass that executable to `AbacusProfile(command=...)`. The command can include an MPI launcher, for example `mpirun -np 4 /path/to/abacus`; ABACUS rank 0 opens the socket connection and broadcasts the i-PI data to the other ranks internally. The ASE interface can be installed from this repository with: ```bash cd interfaces/ASE_interface @@ -153,7 +153,14 @@ with abacus as calc: BFGS(atoms).run(fmax=0.05) ``` -`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. The socket endpoint is selected by the calculator arguments and passed to ABACUS through `ABACUS_SOCKET_ADDRESS`. `unixsocket="abacus_si"` means that ASE listens on the local Unix-domain socket `/tmp/ipi_abacus_si`; abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=/tmp/ipi_abacus_si:UNIX`. For TCP sockets, use the `port` argument instead, which maps to an address such as `localhost:31415`. Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. +`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. It also selects the socket endpoint and passes it to ABACUS through `ABACUS_SOCKET_ADDRESS`, so users normally do not set this environment variable by hand when using abacuslite. + +There are two endpoint styles: + +- `unixsocket="abacus_si"` uses a local Unix-domain socket. ASE creates and listens on `/tmp/ipi_abacus_si`; abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=/tmp/ipi_abacus_si:UNIX`. The `:UNIX` suffix is part of ABACUS' address syntax and means that `/tmp/ipi_abacus_si` is a filesystem socket path, not a TCP host. This is usually the best choice when ASE and ABACUS run on the same node because it avoids TCP port conflicts. +- `port=31415` uses a TCP socket. abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=localhost:31415`, meaning host `localhost` and TCP port `31415`. Use this style when the socket server should listen on a TCP port. If ABACUS is launched manually instead of through `AbacusSocketIO`, set `ABACUS_SOCKET_ADDRESS` yourself to the same `host:port` or `path:UNIX` endpoint. + +Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 74127d1b841..6722db689bb 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -44,7 +44,10 @@ parameters: description: | If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. - [NOTE] Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used. + [NOTE] Use calculation = scf with socket_driver = True. ABACUS connects to the external i-PI server selected by ABACUS_SOCKET_ADDRESS. If ABACUS_SOCKET_ADDRESS is unset, ABACUS uses localhost:31415. The value can use one of two forms: + * host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. + * path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. + When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. default_value: "False" unit: "" availability: "" diff --git a/source/source_io/module_parameter/read_input_item_system.cpp b/source/source_io/module_parameter/read_input_item_system.cpp index bc2d1c94d4f..10444a1b4bd 100644 --- a/source/source_io/module_parameter/read_input_item_system.cpp +++ b/source/source_io/module_parameter/read_input_item_system.cpp @@ -143,7 +143,10 @@ void ReadInput::item_system() item.type = "Boolean"; item.description = R"(If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. -[NOTE] Use calculation = scf with socket_driver = True. Set ABACUS_SOCKET_ADDRESS to host:port or path:UNIX to choose the socket endpoint; if unset, localhost:31415 is used.)"; +[NOTE] Use calculation = scf with socket_driver = True. ABACUS connects to the external i-PI server selected by ABACUS_SOCKET_ADDRESS. If ABACUS_SOCKET_ADDRESS is unset, ABACUS uses localhost:31415. The value can use one of two forms: +* host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. +* path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. +When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument.)"; item.default_value = "False"; read_sync_bool(input.socket_driver); item.check_value = [](const Input_Item& item, const Parameter& para) { From 7c5c0c609f7de64de2481ca061d84272a6ecfd99 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 19 Jul 2026 14:57:32 +0800 Subject: [PATCH 08/24] Fix socket docs generation and launcher parsing --- docs/advanced/input_files/input-main.md | 2 -- docs/advanced/interface/ase.md | 6 +++++- interfaces/ASE_interface/abacuslite/core.py | 21 ++++++++++++++++----- 3 files changed, 21 insertions(+), 8 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index dee59db41d1..a6f71db14cf 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -596,14 +596,12 @@ ### socket_driver - **Type**: Boolean -- **Availability**: *calculation==scf* - **Description**: If set to True, ABACUS keeps the calculation type as scf and receives atomic positions from an external driver through the i-PI socket protocol. > Note: Use calculation = scf with socket_driver = True. ABACUS connects to the external i-PI server selected by ABACUS_SOCKET_ADDRESS. If ABACUS_SOCKET_ADDRESS is unset, ABACUS uses localhost:31415. The value can use one of two forms: - host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. - path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. - When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. - **Default**: False diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index 7e4ed1a6af5..d69fdd33d5f 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -119,7 +119,11 @@ cmake -S . -B build-lcao -DENABLE_MPI=ON -DENABLE_LCAO=ON cmake --build build-lcao --target abacus_basic_para -j ``` -With the ABACUS toolchain workflow, build the normal ABACUS executable with LCAO support when `basis_type=lcao` is needed, then pass that executable to `AbacusProfile(command=...)`. The command can include an MPI launcher, for example `mpirun -np 4 /path/to/abacus`; ABACUS rank 0 opens the socket connection and broadcasts the i-PI data to the other ranks internally. The ASE interface can be installed from this repository with: +With the ABACUS toolchain workflow, build the normal ABACUS executable with LCAO support when `basis_type=lcao` is needed, then pass that executable to `AbacusProfile(command=...)`. The command can include an MPI launcher, for example `mpirun -np 4 /path/to/abacus`; ABACUS rank 0 opens the socket connection and broadcasts the i-PI data to the other ranks internally. On managed clusters, keep scheduler-specific launch options outside the calculator when possible and test the exact launcher command on a compute node. + +For PW calculations on CUDA/ROCm with multiple MPI ranks, use a k-point layout compatible with ABACUS' GPU parallelization. In practice, make sure each k-point pool contains one MPI rank; for example, a 4-rank PW GPU socket calculation should use at least four k-points so the default GPU `kpar` adjustment can assign one rank per pool. A one-k-point PW GPU job with several MPI ranks can fail in the PW GPU transform path; reduce the rank count or use a denser k-point mesh such as a smaller `kspacing`. + +The ASE interface can be installed from this repository with: ```bash cd interfaces/ASE_interface diff --git a/interfaces/ASE_interface/abacuslite/core.py b/interfaces/ASE_interface/abacuslite/core.py index e2c00d4c889..6285b426e5e 100644 --- a/interfaces/ASE_interface/abacuslite/core.py +++ b/interfaces/ASE_interface/abacuslite/core.py @@ -125,11 +125,14 @@ def __init__(self, @staticmethod def parse_version(stdout) -> str: - # up to the ABACUS version v3.9.0.17, the run of command - # `abacus --version` would returns the information organized - # in the following way: - # ABACUS version v3.9.0.17 - return re.match(r'ABACUS version (\S+)', stdout).group(1) + # MPI launchers may add informational lines before ABACUS output. + match = re.search(r'ABACUS version (\S+)', stdout or '') + if match is None: + raise RuntimeError( + 'Could not parse ABACUS version from command output. ' + 'Expected a line like "ABACUS version vX.Y.Z".' + ) + return match.group(1) def get_calculator_command(self, inputfile) -> List[str]: # because ABACUS run in the folder where there are INPUT files, so the @@ -778,6 +781,14 @@ def test_socketio_rejects_cell_changes(self): with self.assertRaisesRegex(PropertyNotImplementedError, 'fixed-cell'): calc._check_fixed_cell(changed) + def test_parse_version_allows_launcher_noise(self): + stdout = 'launcher info\nABACUS version v3.11.0-beta6\n' + self.assertEqual(AbacusProfile.parse_version(stdout), 'v3.11.0-beta6') + + def test_parse_version_rejects_missing_version(self): + with self.assertRaisesRegex(RuntimeError, 'ABACUS version'): + AbacusProfile.parse_version('launcher failed before abacus started') + def test_calculator_results(self): from ase.build.bulk import bulk silicon = bulk('Si', crystalstructure='diamond', a=5.43) From 81a3885eb732c54840634a1177b51b8887573cc7 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Mon, 20 Jul 2026 13:21:21 +0800 Subject: [PATCH 09/24] Document socket running log behavior --- docs/advanced/interface/ase.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index d69fdd33d5f..b9735279d1c 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -168,6 +168,8 @@ Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. +In socket mode, ABACUS keeps one client process alive. All SCF evaluations produced by the same `AbacusSocketIO` instance are appended to the same `OUT.ABACUS/running_scf.log`, because the ABACUS calculation type remains `scf`. The authoritative per-step energy and force results are returned through the i-PI socket to ASE. Use ASE trajectory and optimizer log files, such as `BFGS(atoms, trajectory="opt.traj", logfile="opt.log")`, when each optimizer or MD step should be saved separately. Treat `running_scf.log` mainly as the ABACUS diagnostic log for the socket client, not as one independent FileIO result per structure. + The i-PI protocol does not transmit element symbols. `AbacusSocketIO` therefore sorts the internal socket atoms with the same first-occurrence species grouping used when writing `STRU`, and maps returned forces back to the original ASE `Atoms` order. This avoids silent force/atom mismatches when structures are read from CIF, extxyz, POSCAR, or other formats whose atom order is not already grouped for ABACUS. Users should not manually reorder atoms for socket I/O; pass the physical ASE `Atoms` object directly to the calculator. A complete fixed-cell validation and benchmark example is available in `interfaces/ASE_interface/examples/socketio.py`. From 666f654f797b5efb8e77beed455bbb168df230c4 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Wed, 5 Aug 2026 21:53:12 +0800 Subject: [PATCH 10/24] docs: design fixed-cell i-PI hardening --- ...6-08-05-fixed-cell-ipi-hardening-design.md | 210 ++++++++++++++++++ 1 file changed, 210 insertions(+) create mode 100644 docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md diff --git a/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md b/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md new file mode 100644 index 00000000000..4addb2822d3 --- /dev/null +++ b/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md @@ -0,0 +1,210 @@ +# Fixed-Cell i-PI Hardening Design + +## Objective + +Harden the fixed-cell ABACUS i-PI socket calculator without adding variable-cell +behavior. The work is isolated on `feature/fixed-cell-ipi-hardening`, based +exactly on commit `8b60f83c3e62af75ad57c4c6c61a52a8acdc4d60`. The existing +`feature/fixed-cell-ipi-socket` pull-request branch remains unchanged until the +hardening branch has passed the agreed verification and the user explicitly +authorizes integration. + +The supported scientific scope is fixed-cell SCF evaluation for PW and LCAO on +CPU and GPU. The calculator must return trustworthy energy and forces. It may +return stress only when real ABACUS stress calculation is explicitly enabled; +it must never manufacture a zero stress result and advertise it as computed. + +## Non-Goals + +- Do not update the cell received in `POSDATA`. +- Do not add ASE `UnitCellFilter`, `FrechetCellFilter`, NPT, or barostat support. +- Do not merge or rebase the variable-cell feature into this branch. +- Do not redesign unrelated relaxation, ESolver, MPI, or ASE interfaces. +- Do not install and rebuild an independent ABACUS dependency stack when the + machine's module-provided libraries and existing builds can be reused. + +## Approach + +Implement the hardening directly on the fixed-cell baseline using test-driven +development. The variable-cell branch may be consulted as an already-tested +reference, but changes are introduced as fixed-cell-specific commits with a +failing test preceding each production change. This gives the original PR a +small, reviewable follow-up series and prevents variable-cell control paths from +leaking into the fixed-cell implementation. + +The implementation is divided into four narrow responsibilities: + +1. `IpiSocket` owns raw wire encoding, exact I/O, payload-size arithmetic, and + transport deadlines. +2. A small socket-frame validation unit owns finite-value, shape, cell inverse, + determinant, condition-number, force, energy, stress, and virial contracts. +3. `Socket_Driver` owns the explicit protocol state machine, fixed-cell policy, + ESolver calls, result publication, and MPI failure coordination. +4. `AbacusSocketIO` owns ASE property advertisement, atom-order mapping, process + launch, and client-side validation of returned arrays. + +Dependencies are passed explicitly. No new `GlobalV`, `GlobalC`, or `PARAM` +control dependency is introduced in the socket code. + +## Wire Contract + +The binary protocol uses: + +- exactly 12 bytes for padded ASCII headers; +- `std::int32_t` for replica identifiers, payload lengths, atom counts, and + extra-byte lengths; +- 8-byte IEEE-754 binary64 for cell, inverse cell, positions, energy, forces, + and virial; +- native byte order, matching the official i-PI/NumPy implementation. Mixed- + endian peers are outside the protocol compatibility target and are documented + as unsupported. + +Compile-time assertions enforce integer width, double width, and IEEE-754 +support. All element-count-to-byte-count conversions are checked before +allocation or I/O. The atom count must equal `ucell.nat` before coordinates are +allocated or read. INIT data is capped at 1 MiB because ABACUS does not consume +the initialization payload. + +`IpiSocket` uses bounded connect, receive, and send operations. Timeout values +are supplied explicitly by `Socket_Driver`; they are not mutable workflow state +inside the socket. An operational environment variable, +`ABACUS_SOCKET_TIMEOUT_SECONDS`, selects a positive finite timeout and defaults +to 300 seconds. It applies only while waiting for socket transport, not while +ABACUS is running SCF. The existing documentation is updated with the default, +accepted values, and failure behavior. + +## Protocol State Machine + +Use an explicit local enum with three states: + +```text +NeedInit --INIT--> Ready --POSDATA/compute--> HaveData --GETFORCE--> Ready +``` + +`STATUS` is valid in every state and returns `NEEDINIT`, `READY`, or `HAVEDATA`. +`EXIT` is accepted in `NeedInit` or `Ready`. A clean peer close is accepted only +when no computed frame is pending. The following are protocol errors: + +- INIT outside `NeedInit`; +- POSDATA outside `Ready`; +- GETFORCE outside `HaveData`; +- EXIT or peer close while a computed frame is pending; +- unknown or malformed headers. + +No default, stale, or partially computed frame is published after an error. + +## Numerical and Scientific Contract + +Before changing ABACUS state, validate the complete input frame: + +- all cell, inverse-cell, and position entries are finite binary64 values; +- the cell determinant is finite and positive; +- the cell condition number is below `1.0e12`; +- the received inverse is consistent with the received cell using combined + absolute and condition-scaled tolerances; +- the atom count equals the STRU atom count; +- the received fixed cell matches the initial ABACUS cell using a combined + absolute/relative tolerance that remains strict at normal cell sizes. + +The wire matrix is the i-PI column-vector matrix `H`. ASE sends its row-vector +matrix `A` as `A.T`, while official i-PI sends `H` directly. The same physical +cell therefore obeys `A = H.T`. Tests use a nonsymmetric triclinic cell so a +mistaken transpose cannot pass accidentally. + +After `runner()`: + +- reject a KS frame when SCF did not converge; +- require finite energy; +- require force shape `nat x 3` and finite force entries; +- when stress is enabled, require a finite `3 x 3` stress tensor, positive + finite volume, and a sufficiently symmetric stress tensor; +- convert ABACUS Ry-based energy, force, and stress to i-PI Hartree/Bohr wire + units with explicit constants; +- construct wire virial with the verified sign and transpose convention; +- mark `HaveData` only after the full output frame passes validation. + +The initial POSDATA atom order must be PBC-equivalent to STRU order. Because the +wire protocol carries no species, a mismatch is a fatal contract violation for +raw i-PI clients rather than a warning. The ASE wrapper continues to sort atoms +into STRU species order and maps returned forces back to the caller's order. + +## ASE Property Semantics + +`AbacusSocketIO` defines instance-level `implemented_properties`: + +- default fixed-cell mode: `energy`, `free_energy`, and `forces`; +- fixed-cell mode with `inp['cal_stress'] = 1`: the same properties plus + `stress`. + +When stress is enabled, `_socket_inp` preserves/enforces `cal_stress=1`, the C++ +driver computes real ABACUS stress, and the wrapper converts the returned virial +to ASE stress. When it is disabled, asking ASE for stress raises +`PropertyNotImplementedError`; zero virial is not exposed as a physical result. + +The wrapper validates result energy, forces, and virial for shape and finiteness +before updating `self.results`. Python `assert` is not relied upon for runtime +validation because optimized Python may disable assertions. + +## MPI Failure Handling + +Only rank zero performs socket I/O. Root I/O failures are broadcast before any +rank advances. For each collective computation stage—position application, +ESolver execution, energy/force/stress extraction, and result validation—local +exceptions are collected across ranks before the next collective operation. + +If one rank fails while others may be inside or approaching an incompatible +collective, terminate the socket calculation consistently using the established +ABACUS MPI failure mechanism. Focused MPI tests inject failures on root and +non-root ranks and require bounded termination rather than a hang. + +## Test Strategy + +Every behavior change follows red-green-refactor and is committed separately. +The required automated coverage is: + +1. Wire unit tests: exact 12-byte headers, exact four-byte integers, exact + eight-byte doubles, split TCP/UNIX payloads, clean and partial closes, + overflow rejection, INIT-size rejection, and transport timeout. +2. Frame unit tests: NaN/Inf, singular/negative/ill-conditioned cells, + inconsistent inverse, atom-count overflow/mismatch, changed fixed cell, + invalid energy/force/stress, stress symmetry, unit conversion, virial sign, + and transpose convention. +3. Driver tests: complete valid sequence and every invalid state transition, + nonconverged SCF, child cleanup, unknown/EXIT headers, pending-frame close, + and no publication of invalid or stale frames. +4. MPI tests: root and non-root injected failures terminate within a fixed + deadline and report the failing stage. +5. ASE tests: property advertisement, stress gating, finite/shape validation, + species sorting, force remapping, fixed-cell rejection, and process cleanup. +6. Golden-wire tests: ASE and official i-PI 3.2.0 produce the same wire `H` for + the same nonsymmetric physical cell `A = H.T`. +7. Runtime calculations, with `OMP_NUM_THREADS=1`: one-frame and short repeated- + frame comparisons against ABACUS FileIO for CPU PW and CPU LCAO, followed by + one-GPU PW and one-GPU LCAO checks. Energy, forces, and optional stress use + explicit tolerances recorded with each result. + +Runtime builds reuse the locally available ABACUS development modules and linked +libraries. Only the changed targets are rebuilt where CMake permits it. MPI and +socket runtime tests run outside the restricted sandbox. + +## Verification Gates + +The hardening branch is ready to propose for integration only when all of the +following are true: + +- focused C++/Python/MPI tests pass from a fresh build of the changed targets; +- existing socket and INPUT parser tests pass; +- CPU PW and LCAO energy/force comparisons pass; +- official i-PI 3.2.0 and ASE nonsymmetric-cell contract tests pass; +- one-GPU PW and LCAO socket calculations pass, or an external hardware blocker + is reported without claiming completion; +- stress is either verified as real with sign/unit/format checks or is correctly + unavailable; +- `git diff --check`, ABACUS agent governance checks, and relevant pre-commit + checks pass or have concrete review rationales; +- the worktree is clean after committed changes; +- the original PR branch SHA remains unchanged. + +After these gates, push only `feature/fixed-cell-ipi-hardening` and report its +commit SHA and exact verification evidence. Merging or pushing into +`feature/fixed-cell-ipi-socket` requires a separate explicit user authorization. From bfc028bbf19d0253c195249760717fb69b7f4760 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Wed, 12 Aug 2026 19:25:49 +0800 Subject: [PATCH 11/24] docs: record ELPA 2026.02 genELPA compatibility issue --- .../elpa_2026_02_genelpa_issue.md | 190 ++++++++++++++++++ docs/developers_guide/index.rst | 1 + 2 files changed, 191 insertions(+) create mode 100644 docs/developers_guide/elpa_2026_02_genelpa_issue.md diff --git a/docs/developers_guide/elpa_2026_02_genelpa_issue.md b/docs/developers_guide/elpa_2026_02_genelpa_issue.md new file mode 100644 index 00000000000..56571cccbf0 --- /dev/null +++ b/docs/developers_guide/elpa_2026_02_genelpa_issue.md @@ -0,0 +1,190 @@ +# ELPA 2026.02.001 `genelpa` MPI/grid compatibility report + +This developer note records a reproducible external-solver issue found while +validating the ABACUS i-PI socket calculator. It is separate from the socket +implementation PRs: the socket code does not cause this failure. + +## Summary + +The installed ELPA 2026.02.001 build reproducibly fails for a minimal public +ELPA C API program using a complex-double, two-stage eigensolver on a row-major +`1 x 2` BLACS process grid. The same program passes with ELPA 2026.02.001 on a +`2 x 1` grid and with ELPA 2025.06.001 on both grids. The one-stage solver +passes in all tested ELPA 2026.02.001 cases. + +The evidence is consistent with an ELPA 2026.02.001 regression or +build-specific incompatibility in this multi-process/grid configuration. It is +not yet sufficient to claim that an upstream ELPA bug has been confirmed; the +reproducer should be checked against another 2026.02.001 build or by ELPA +maintainers. + +## ABACUS configuration + +The failure was first observed in an ABACUS LCAO calculation with: + +```text +nlocal = 26 +nbands = 14 +nb2d = 1 +MPI ranks = 2 +row-major BLACS grid = 1 x 2 +complex Hamiltonian path +ks_solver = genelpa +``` + +For two MPI ranks, ABACUS's row-major `Parallel_2D` setup produces the `1 x 2` +grid. For matrices smaller than 500 orbitals, the default LCAO block size is +`nb2d = 1`. + +Relevant code paths are: + +```text +source/source_lcao/LCAO_init_basis.cpp +source/source_base/parallel_2d.cpp +source/source_hsolver/module_genelpa/elpa_new.cpp +source/source_hsolver/module_genelpa/elpa_new_complex.cpp +``` + +## Standalone reproducer + +The reproducer calls the ELPA C API directly and does not call ABACUS or the +i-PI socket layer. It uses a small positive-definite Hermitian complex matrix: + +```text +global matrix size: n = 26 +number of eigenvectors: nev = 14 +block size: nblk = 1 +MPI ranks: 2 +BLACS layout: row-major +process grids: 1 x 2 and 2 x 1 +``` + +The temporary reproducer was `/tmp/elpa_order_test.cpp`. The essential command +was: + +```bash +OMP_NUM_THREADS=1 \ +OMPI_MCA_pml=ob1 \ +OMPI_MCA_btl=self,tcp \ +mpirun --oversubscribe -np 2 ./elpa_order_test \ + two_before b1 1x2 26 1 +``` + +The test was repeated with solver options set before and after `elpa_setup()`, +with ELPA 2026.02.001 and 2025.06.001, and with several complex kernels. + +## Result matrix + +`RC=0` means that the program completed and ELPA returned `ELPA_OK` on all +ranks. Non-zero results are launcher or process failures. + +| ELPA version | BLACS grid | one-stage | two-stage | +| --- | ---: | ---: | ---: | +| 2026.02.001 | `1 x 2` | pass (`RC=0`) | fail | +| 2026.02.001 | `2 x 1` | pass (`RC=0`) | pass (`RC=0`) | +| 2025.06.001 | `1 x 2` | pass (`RC=0`) | pass (`RC=0`) | +| 2025.06.001 | `2 x 1` | pass (`RC=0`) | pass (`RC=0`) | + +The ELPA 2026.02.001 `1 x 2` two-stage failure was reproduced with: + +```text +ELPA_2STAGE_COMPLEX_AVX512_BLOCK1 +ELPA_2STAGE_COMPLEX_AVX512_BLOCK2 +ELPA_2STAGE_COMPLEX_GENERIC +ELPA_2STAGE_COMPLEX_GENERIC_SIMPLE +``` + +Setting `solver` and `complex_kernel` before or after `elpa_setup()` did not +change the result. Linking the same reproducer with OpenMPI 5.0.8 instead of +OpenMPI 5.0.10 also reproduced the failure. + +## Observed errors + +With OpenMPI OB1/TCP, representative output is: + +```text +mca_btl_tcp_frag_send: writev error +Bad address(3) +An error occurred in Socket closed +MPI_ERRORS_ARE_FATAL +``` + +With default UCX, representative runs reported UCX/xpmem errors followed by +`SIGBUS` or `SIGSEGV`. The stack reached: + +```text +MPI_Bcast + -> trans_ev_tridi_to_band_complex_double + -> elpa_solve_evp_complex_2stage + -> elpa_eigenvectors +``` + +Some runs also reported: + +```text +Operating system error: Cannot allocate memory +Integer overflow in xmallocarray +``` + +These are MPI/ELPA internal errors, not i-PI protocol socket errors. + +## Current attribution + +The experiments support these conclusions: + +1. The issue is reproducible without ABACUS and without the i-PI socket layer. +2. It is not explained by the `gint_precision` input keyword. That keyword is + unsupported by the LTS 3.10.1 parser, but removing it does not fix the + develop-module failure. +3. It is not only an UCX/xpmem problem: OB1/TCP runs fail as well. +4. It is not caused by one selected complex kernel or by setting the solver + option before versus after `elpa_setup()`. +5. The version/grid/algorithm matrix points to an ELPA 2026.02.001 two-stage, + multi-process compatibility problem. + +The safe issue wording is: + +> ELPA 2026.02.001 reproducibly fails for complex-double two-stage +> diagonalization on an ABACUS-compatible `1 x 2` BLACS grid; ELPA 2025.06.001 +> passes the same reproducer. + +Do not state that an upstream ELPA bug is confirmed until another 2026.02.001 +build or ELPA maintainers reproduce the issue. + +## ABACUS wrapper observations + +The ABACUS `module_genelpa` wrapper has separate engineering issues that +should be handled in a solver-focused change, not hidden inside the socket +calculator PR: + +- `elpa_set()` and `elpa_setup()` return codes are not consistently checked; +- some `MPI_Allreduce()` results are subsequently overwritten by local `info`; +- ELPA handle and `elpa_init()`/`elpa_uninit()` lifetimes need a focused audit; +- failures should become rank-consistent, actionable ABACUS diagnostics. + +These observations improve ABACUS diagnostics, but do not explain why a +standalone public-API program reproduces the ELPA failure. + +## Temporary workarounds + +Until the dependency issue is resolved, use one of the following for affected +CPU LCAO runs: + +```text +ELPA 2025.06.001 +ks_solver = elpa +ks_solver = scalapack_gvx +``` + +Socket validation should use a solver configuration known to work on the +target machine. Socket PRs should not silently override ELPA behavior. + +## Follow-up + +1. Submit the reproducer and result matrix to ELPA maintainers. +2. Repeat with an independent ELPA 2026.02.001 CPU-only build and on a + GPU-capable compute node. +3. If confirmed externally, document the affected ELPA/build/grid combination + as a platform dependency limitation. +4. Independently harden the ABACUS genELPA wrapper with checked return codes + and collective error propagation. diff --git a/docs/developers_guide/index.rst b/docs/developers_guide/index.rst index f346894be47..1807f2dc77d 100644 --- a/docs/developers_guide/index.rst +++ b/docs/developers_guide/index.rst @@ -11,3 +11,4 @@ This section provides guidelines and resources for developers working on the ABA basic_types_class.md agent_governance.md + elpa_2026_02_genelpa_issue.md From 4dc207e19e98fa268f323f90f89dc60f69257599 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Tue, 4 Aug 2026 00:28:59 +0800 Subject: [PATCH 12/24] fix: use explicit i-PI wire integer type --- source/source_relax/socket_driver.cpp | 10 +- source/source_relax/socket_ipi.cpp | 27 ++- source/source_relax/socket_ipi.h | 5 +- source/source_relax/test/socket_ipi_test.cpp | 189 +++++++++++++++++++ 4 files changed, 219 insertions(+), 12 deletions(-) diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index 42fafda5e63..b97e4c118d2 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -371,8 +371,8 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { try { - rid = socket.read_int(); - nbytes = socket.read_int(); + rid = socket.read_int32(); + nbytes = socket.read_int32(); if (nbytes < 0) { io_failed = 1; @@ -416,7 +416,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { cell = socket.read_doubles(9); inv_cell = socket.read_doubles(9); - nat_socket = socket.read_int(); + nat_socket = socket.read_int32(); if (nat_socket < 0) { io_failed = 1; @@ -494,10 +494,10 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { socket.write_header("FORCEREADY"); socket.write_double(energy_hartree); - socket.write_int(nat_return); + socket.write_int32(nat_return); socket.write_doubles(forces_hartree_bohr); socket.write_doubles(virial_hartree); - socket.write_int(0); + socket.write_int32(0); } catch (const std::exception& exc) { diff --git a/source/source_relax/socket_ipi.cpp b/source/source_relax/socket_ipi.cpp index 2bb09385b0d..1c2788dae6c 100644 --- a/source/source_relax/socket_ipi.cpp +++ b/source/source_relax/socket_ipi.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -10,6 +11,11 @@ #include #include +static_assert(sizeof(std::int32_t) == 4, "i-PI requires a 4-byte integer"); +static_assert(sizeof(double) == 8, "i-PI requires an 8-byte float"); +static_assert(std::numeric_limits::is_iec559, + "i-PI requires IEEE-754 double precision"); + namespace { constexpr std::size_t IPI_HEADER_LEN = 12; @@ -39,6 +45,15 @@ std::string padded_header(const std::string& header) out.resize(IPI_HEADER_LEN, ' '); return out; } + +std::size_t checked_double_bytes(std::size_t n) +{ + if (n > SIZE_MAX / sizeof(double)) + { + throw std::overflow_error("i-PI double payload byte count overflows for " + std::to_string(n) + " elements"); + } + return n * sizeof(double); +} } // namespace IpiSocketClosed::IpiSocketClosed(const std::string& message) : std::runtime_error(message) @@ -166,14 +181,14 @@ void IpiSocket::write_header(const std::string& header) this->write_exact(padded.data(), padded.size()); } -int IpiSocket::read_int() +std::int32_t IpiSocket::read_int32() { - int value = 0; + std::int32_t value = 0; this->read_exact(&value, sizeof(value)); return value; } -void IpiSocket::write_int(int value) +void IpiSocket::write_int32(std::int32_t value) { this->write_exact(&value, sizeof(value)); } @@ -192,19 +207,21 @@ void IpiSocket::write_double(double value) std::vector IpiSocket::read_doubles(std::size_t n) { + const std::size_t nbytes = checked_double_bytes(n); std::vector values(n); if (!values.empty()) { - this->read_exact(values.data(), values.size() * sizeof(double)); + this->read_exact(values.data(), nbytes); } return values; } void IpiSocket::write_doubles(const std::vector& values) { + const std::size_t nbytes = checked_double_bytes(values.size()); if (!values.empty()) { - this->write_exact(values.data(), values.size() * sizeof(double)); + this->write_exact(values.data(), nbytes); } } diff --git a/source/source_relax/socket_ipi.h b/source/source_relax/socket_ipi.h index eac15eb41ac..563a111f95c 100644 --- a/source/source_relax/socket_ipi.h +++ b/source/source_relax/socket_ipi.h @@ -2,6 +2,7 @@ #define ABACUS_SOCKET_IPI_H #include +#include #include #include #include @@ -27,8 +28,8 @@ class IpiSocket std::string read_header(); void write_header(const std::string& header); - int read_int(); - void write_int(int value); + std::int32_t read_int32(); + void write_int32(std::int32_t value); double read_double(); void write_double(double value); diff --git a/source/source_relax/test/socket_ipi_test.cpp b/source/source_relax/test/socket_ipi_test.cpp index 930ab1fe629..cbe92707560 100644 --- a/source/source_relax/test/socket_ipi_test.cpp +++ b/source/source_relax/test/socket_ipi_test.cpp @@ -3,15 +3,18 @@ #include "gtest/gtest.h" #include +#include #include #include #include +#include #include #include #include #include #include #include +#include namespace { @@ -254,3 +257,189 @@ TEST(IpiSocketTest, PartialHeaderCloseStaysRuntimeError) peer.join(); rethrow_thread_error(thread_error); } + +TEST(IpiSocketTest, Int32UsesExactlyFourNativeEndianBytes) +{ + UnixSocketServer server; + const std::int32_t expected = INT32_C(0x12345678); + std::vector received(4); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + recv_all(fd, received.data(), received.size()); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + socket.write_int32(expected); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); + EXPECT_EQ(0, std::memcmp(received.data(), &expected, 4)); +} + +TEST(IpiSocketTest, DoubleUsesExactlyEightNativeEndianBytes) +{ + UnixSocketServer server; + const double expected = -1234.5; + std::vector received(8); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + recv_all(fd, received.data(), received.size()); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + socket.write_double(expected); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); + EXPECT_EQ(0, std::memcmp(received.data(), &expected, 8)); +} + +TEST(IpiSocketTest, ReadInt32HandlesSplitPayload) +{ + UnixSocketServer server; + const std::int32_t expected = INT32_C(0x12345678); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + const char* bytes = reinterpret_cast(&expected); + send_all(fd, bytes, 2); + send_all(fd, bytes + 2, 2); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + EXPECT_EQ(expected, socket.read_int32()); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); +} + +TEST(IpiSocketTest, ReadInt32RejectsMidPayloadClose) +{ + UnixSocketServer server; + const std::int32_t value = INT32_C(0x12345678); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + send_all(fd, &value, 2); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + EXPECT_THROW(socket.read_int32(), IpiSocketClosed); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); +} + +TEST(IpiSocketTest, WriteDoublesCompletesLargePayloadWithSmallPeerReads) +{ + UnixSocketServer server; + std::vector expected(1 << 18); + for (std::size_t i = 0; i < expected.size(); ++i) + { + expected[i] = -1234.5 + static_cast(i) * 0.25; + } + std::vector received(expected.size() * sizeof(double)); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + std::size_t done = 0; + while (done < received.size()) + { + const std::size_t remaining = received.size() - done; + const std::size_t chunk = remaining < 37 ? remaining : 37; + const ssize_t nread = ::recv(fd, received.data() + done, chunk, 0); + if (nread < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("recv failed")); + } + if (nread == 0) + { + throw std::runtime_error("socket closed while receiving large payload"); + } + done += static_cast(nread); + } + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + socket.write_doubles(expected); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); + EXPECT_EQ(0, std::memcmp(received.data(), expected.data(), received.size())); +} + +TEST(IpiSocketTest, ReadDoublesRejectsByteCountOverflow) +{ + IpiSocket socket; + const std::size_t count = std::numeric_limits::max() / sizeof(double) + 1; + + try + { + static_cast(socket.read_doubles(count)); + FAIL() << "overflowing double payload size should throw"; + } + catch (const std::overflow_error& exc) + { + EXPECT_NE(std::string::npos, std::string(exc.what()).find(std::to_string(count))); + } + catch (...) + { + FAIL() << "overflowing double payload size should throw std::overflow_error"; + } +} From 9b63794324efbadafbf6e127060bb8a505c19842 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Tue, 4 Aug 2026 00:45:55 +0800 Subject: [PATCH 13/24] feat: add socket frame numerical contract --- source/source_relax/CMakeLists.txt | 1 + source/source_relax/socket_frame.cpp | 418 ++++++++++++++++++ source/source_relax/socket_frame.h | 51 +++ source/source_relax/test/CMakeLists.txt | 5 + .../source_relax/test/socket_frame_test.cpp | 278 ++++++++++++ 5 files changed, 753 insertions(+) create mode 100644 source/source_relax/socket_frame.cpp create mode 100644 source/source_relax/socket_frame.h create mode 100644 source/source_relax/test/socket_frame_test.cpp diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 77079eb6a1f..f366af77beb 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -3,6 +3,7 @@ add_library( OBJECT relax_data.cpp socket_ipi.cpp + socket_frame.cpp socket_driver.cpp cg_base.cpp relax_driver.cpp diff --git a/source/source_relax/socket_frame.cpp b/source/source_relax/socket_frame.cpp new file mode 100644 index 00000000000..466499c198e --- /dev/null +++ b/source/source_relax/socket_frame.cpp @@ -0,0 +1,418 @@ +#include "socket_frame.h" + +#include +#include +#include + +namespace +{ +const int MATRIX_DIMENSION = 3; +const int MAX_JACOBI_SWEEPS = 32; + +bool is_finite_matrix(const SocketFrame::Matrix9& values) +{ + for (std::size_t index = 0; index < values.size(); ++index) + { + if (!std::isfinite(values[index])) + { + return false; + } + } + return true; +} + +double column_norm_squared(const SocketFrame::Matrix9& values, int column) +{ + double norm_squared = 0.0; + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + const double value = values[row * MATRIX_DIMENSION + column]; + norm_squared += value * value; + } + return norm_squared; +} + +double column_dot(const SocketFrame::Matrix9& values, int first, int second) +{ + double dot = 0.0; + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + dot += values[row * MATRIX_DIMENSION + first] * values[row * MATRIX_DIMENSION + second]; + } + return dot; +} + +bool columns_are_orthogonal(const SocketFrame::Matrix9& values) +{ + const double multiplier = 32.0 * std::numeric_limits::epsilon(); + const int pairs[3][2] = {{0, 1}, {0, 2}, {1, 2}}; + for (int pair = 0; pair < 3; ++pair) + { + const int first = pairs[pair][0]; + const int second = pairs[pair][1]; + const double first_norm = column_norm_squared(values, first); + const double second_norm = column_norm_squared(values, second); + const double tolerance = multiplier * std::sqrt(first_norm * second_norm); + if (std::fabs(column_dot(values, first, second)) > tolerance) + { + return false; + } + } + return true; +} + +void rotate_columns(SocketFrame::Matrix9& values, int first, int second, double cosine, double sine) +{ + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + const int first_index = row * MATRIX_DIMENSION + first; + const int second_index = row * MATRIX_DIMENSION + second; + const double first_value = values[first_index]; + const double second_value = values[second_index]; + values[first_index] = cosine * first_value - sine * second_value; + values[second_index] = sine * first_value + cosine * second_value; + } +} + +bool one_sided_jacobi(SocketFrame::Matrix9& columns, SocketFrame::Matrix9& right_vectors) +{ + right_vectors = {{1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0}}; + const double multiplier = 32.0 * std::numeric_limits::epsilon(); + const int pairs[3][2] = {{0, 1}, {0, 2}, {1, 2}}; + + for (int sweep = 0; sweep < MAX_JACOBI_SWEEPS; ++sweep) + { + for (int pair = 0; pair < 3; ++pair) + { + const int first = pairs[pair][0]; + const int second = pairs[pair][1]; + const double first_norm = column_norm_squared(columns, first); + const double second_norm = column_norm_squared(columns, second); + const double dot = column_dot(columns, first, second); + const double tolerance = multiplier * std::sqrt(first_norm * second_norm); + if (std::fabs(dot) <= tolerance) + { + continue; + } + + const double tau = (second_norm - first_norm) / (2.0 * dot); + const double tangent + = std::copysign(1.0 / (std::fabs(tau) + std::hypot(1.0, tau)), tau); + const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent); + const double sine = tangent * cosine; + rotate_columns(columns, first, second, cosine, sine); + rotate_columns(right_vectors, first, second, cosine, sine); + } + + if (columns_are_orthogonal(columns)) + { + return true; + } + } + return false; +} + +long double scaled_determinant(const SocketFrame::Matrix9& values) +{ + const long double a00 = values[0]; + const long double a01 = values[1]; + const long double a02 = values[2]; + const long double a10 = values[3]; + const long double a11 = values[4]; + const long double a12 = values[5]; + const long double a20 = values[6]; + const long double a21 = values[7]; + const long double a22 = values[8]; + return a00 * (a11 * a22 - a12 * a21) + - a01 * (a10 * a22 - a12 * a20) + + a02 * (a10 * a21 - a11 * a20); +} + +double received_inverse_residual(const SocketFrame::Matrix9& cell, + const SocketFrame::Matrix9& inverse) +{ + long double maximum = 0.0L; + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + for (int column = 0; column < MATRIX_DIMENSION; ++column) + { + long double product = 0.0L; + for (int inner = 0; inner < MATRIX_DIMENSION; ++inner) + { + product += static_cast(cell[row * MATRIX_DIMENSION + inner]) + * inverse[inner * MATRIX_DIMENSION + column]; + } + const long double expected = row == column ? 1.0L : 0.0L; + maximum = std::max(maximum, std::fabs(product - expected)); + } + } + return static_cast(maximum); +} +} // namespace + +namespace SocketFrame +{ +Matrix9 transpose_matrix9(const Matrix9& values) +{ + return {{values[0], values[3], values[6], + values[1], values[4], values[7], + values[2], values[5], values[8]}}; +} + +CellValidation validate_ipi_cell(const Matrix9& cell_wire, + const Matrix9& inverse_wire, + double max_condition_number, + double inverse_absolute_tolerance, + double inverse_relative_tolerance) +{ + CellValidation result; + result.ok = false; + result.message.clear(); + result.determinant_bohr3 = 0.0; + result.condition_number_2 = std::numeric_limits::infinity(); + result.inverse_residual = std::numeric_limits::infinity(); + result.computed_inverse_wire_bohr_inv.fill(0.0); + + if (!is_finite_matrix(cell_wire) || !is_finite_matrix(inverse_wire)) + { + result.message = "cell and received inverse entries must be finite"; + return result; + } + if (!std::isfinite(max_condition_number) || max_condition_number <= 0.0 + || !std::isfinite(inverse_absolute_tolerance) || inverse_absolute_tolerance < 0.0 + || !std::isfinite(inverse_relative_tolerance) || inverse_relative_tolerance < 0.0) + { + result.message = "cell validation tolerances must be finite and nonnegative"; + return result; + } + + double scale = 0.0; + for (std::size_t index = 0; index < cell_wire.size(); ++index) + { + scale = std::max(scale, std::fabs(cell_wire[index])); + } + if (scale == 0.0) + { + result.message = "cell determinant must be positive"; + return result; + } + + Matrix9 scaled_cell; + for (std::size_t index = 0; index < cell_wire.size(); ++index) + { + scaled_cell[index] = cell_wire[index] / scale; + } + const long double determinant_scaled = scaled_determinant(scaled_cell); + if (determinant_scaled <= 0.0L) + { + result.message = "cell determinant must be positive"; + return result; + } + const long double scale_long = scale; + const long double determinant + = determinant_scaled * scale_long * scale_long * scale_long; + if (!std::isfinite(determinant) + || determinant > static_cast(std::numeric_limits::max())) + { + result.message = "cell determinant is not representable as a finite double"; + return result; + } + result.determinant_bohr3 = static_cast(determinant); + if (!std::isfinite(result.determinant_bohr3) || result.determinant_bohr3 <= 0.0) + { + result.message = "cell determinant is not representable as a positive finite double"; + return result; + } + + Matrix9 orthogonal_columns = scaled_cell; + Matrix9 right_vectors; + if (!one_sided_jacobi(orthogonal_columns, right_vectors)) + { + result.message = "cell singular-value iteration did not converge"; + return result; + } + + double singular_values[MATRIX_DIMENSION]; + double largest_singular = 0.0; + double smallest_singular = std::numeric_limits::infinity(); + for (int column = 0; column < MATRIX_DIMENSION; ++column) + { + singular_values[column] = std::sqrt(column_norm_squared(orthogonal_columns, column)); + largest_singular = std::max(largest_singular, singular_values[column]); + smallest_singular = std::min(smallest_singular, singular_values[column]); + } + if (smallest_singular == 0.0 || !std::isfinite(smallest_singular)) + { + result.message = "cell is singular"; + return result; + } + result.condition_number_2 = largest_singular / smallest_singular; + if (!std::isfinite(result.condition_number_2) + || result.condition_number_2 > max_condition_number) + { + result.message = "cell condition number exceeds the configured maximum"; + return result; + } + + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + for (int column = 0; column < MATRIX_DIMENSION; ++column) + { + long double inverse_value = 0.0L; + for (int singular = 0; singular < MATRIX_DIMENSION; ++singular) + { + const long double sigma = singular_values[singular]; + inverse_value + += static_cast(right_vectors[row * MATRIX_DIMENSION + singular]) + * orthogonal_columns[column * MATRIX_DIMENSION + singular] + / (static_cast(scale) * sigma * sigma); + } + result.computed_inverse_wire_bohr_inv[row * MATRIX_DIMENSION + column] + = static_cast(inverse_value); + } + } + + result.inverse_residual = received_inverse_residual(cell_wire, inverse_wire); + const double residual_limit + = inverse_absolute_tolerance + + inverse_relative_tolerance * result.condition_number_2 + * std::numeric_limits::epsilon(); + if (!std::isfinite(result.inverse_residual) || result.inverse_residual > residual_limit) + { + result.message = "received cell inverse is inconsistent with the cell"; + return result; + } + + result.ok = true; + return result; +} + +bool validate_positions(const std::vector& positions_bohr, + std::size_t coordinate_count, + std::string& message) +{ + if (positions_bohr.size() != coordinate_count) + { + message = "position coordinate count does not match the validated atom count"; + return false; + } + for (std::size_t index = 0; index < positions_bohr.size(); ++index) + { + if (!std::isfinite(positions_bohr[index])) + { + message = "position coordinates must be finite"; + return false; + } + } + message.clear(); + return true; +} + +bool checked_position_count(std::int32_t nat_socket, + int nat_expected, + std::size_t& coordinate_count, + std::string& message) +{ + if (nat_socket != nat_expected) + { + message = "socket atom count does not match the expected atom count"; + return false; + } + if (nat_socket < 0) + { + message = "socket atom count must not be negative"; + return false; + } + const std::size_t atom_count = static_cast(nat_socket); + if (atom_count > std::numeric_limits::max() / 3) + { + message = "socket position coordinate count is not representable"; + return false; + } + coordinate_count = 3 * atom_count; + message.clear(); + return true; +} + +VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, + double volume_bohr3, + double antisymmetric_absolute_tolerance, + double antisymmetric_relative_tolerance) +{ + VirialConversion result; + result.ok = false; + result.message.clear(); + result.wire_virial_hartree.fill(0.0); + result.max_antisymmetric_component = 0.0; + + if (!is_finite_matrix(stress_ry_per_bohr3)) + { + result.message = "stress entries must be finite"; + return result; + } + if (!std::isfinite(volume_bohr3) || volume_bohr3 <= 0.0) + { + result.message = "cell volume must be finite and positive"; + return result; + } + if (!std::isfinite(antisymmetric_absolute_tolerance) + || antisymmetric_absolute_tolerance < 0.0 + || !std::isfinite(antisymmetric_relative_tolerance) + || antisymmetric_relative_tolerance < 0.0) + { + result.message = "stress symmetry tolerances must be finite and nonnegative"; + return result; + } + + double maximum_stress = 0.0; + for (std::size_t index = 0; index < stress_ry_per_bohr3.size(); ++index) + { + maximum_stress = std::max(maximum_stress, std::fabs(stress_ry_per_bohr3[index])); + } + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + for (int column = row + 1; column < MATRIX_DIMENSION; ++column) + { + const double difference + = std::fabs(stress_ry_per_bohr3[row * MATRIX_DIMENSION + column] + - stress_ry_per_bohr3[column * MATRIX_DIMENSION + row]); + result.max_antisymmetric_component + = std::max(result.max_antisymmetric_component, difference); + } + } + const double symmetry_limit + = antisymmetric_absolute_tolerance + antisymmetric_relative_tolerance * maximum_stress; + if (!std::isfinite(result.max_antisymmetric_component) + || result.max_antisymmetric_component > symmetry_limit) + { + result.message = "stress tensor is not symmetric within tolerance"; + return result; + } + + Matrix9 virial; + for (int row = 0; row < MATRIX_DIMENSION; ++row) + { + for (int column = 0; column < MATRIX_DIMENSION; ++column) + { + const long double symmetric_stress + = 0.5L + * (static_cast(stress_ry_per_bohr3[row * MATRIX_DIMENSION + column]) + + stress_ry_per_bohr3[column * MATRIX_DIMENSION + row]); + const long double converted = 0.5L * volume_bohr3 * symmetric_stress; + if (!std::isfinite(converted) + || std::fabs(converted) + > static_cast(std::numeric_limits::max())) + { + result.message = "converted virial is not representable as finite doubles"; + return result; + } + virial[row * MATRIX_DIMENSION + column] = static_cast(converted); + } + } + result.wire_virial_hartree = transpose_matrix9(virial); + result.ok = true; + return result; +} +} // namespace SocketFrame diff --git a/source/source_relax/socket_frame.h b/source/source_relax/socket_frame.h new file mode 100644 index 00000000000..759a4ee13a2 --- /dev/null +++ b/source/source_relax/socket_frame.h @@ -0,0 +1,51 @@ +#ifndef SOURCE_RELAX_SOCKET_FRAME_H +#define SOURCE_RELAX_SOCKET_FRAME_H + +#include +#include +#include +#include +#include + +namespace SocketFrame +{ +using Matrix9 = std::array; + +struct CellValidation +{ + bool ok; + std::string message; + double determinant_bohr3; + double condition_number_2; + double inverse_residual; + Matrix9 computed_inverse_wire_bohr_inv; +}; + +struct VirialConversion +{ + bool ok; + std::string message; + Matrix9 wire_virial_hartree; + double max_antisymmetric_component; +}; + +Matrix9 transpose_matrix9(const Matrix9& values); +CellValidation validate_ipi_cell(const Matrix9& cell_wire, + const Matrix9& inverse_wire, + double max_condition_number, + double inverse_absolute_tolerance, + double inverse_relative_tolerance); +bool validate_positions(const std::vector& positions_bohr, + std::size_t coordinate_count, + std::string& message); +bool checked_position_count(std::int32_t nat_socket, + int nat_expected, + std::size_t& coordinate_count, + std::string& message); +VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, + double volume_bohr3, + double antisymmetric_absolute_tolerance, + double antisymmetric_relative_tolerance); +} // namespace SocketFrame + +#endif diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 322751c352f..3b50983345a 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -12,6 +12,11 @@ AddTest( SOURCES socket_ipi_test.cpp ../socket_ipi.cpp ) +AddTest( + TARGET MODULE_RELAX_socket_frame_test + SOURCES socket_frame_test.cpp ../socket_frame.cpp +) + AddTest( TARGET MODULE_RELAX_relax_new_line_search LIBS parameter diff --git a/source/source_relax/test/socket_frame_test.cpp b/source/source_relax/test/socket_frame_test.cpp new file mode 100644 index 00000000000..26e25938e58 --- /dev/null +++ b/source/source_relax/test/socket_frame_test.cpp @@ -0,0 +1,278 @@ +#include "../socket_frame.h" + +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include + +namespace +{ +using SocketFrame::CellValidation; +using SocketFrame::Matrix9; +using SocketFrame::VirialConversion; +using SocketFrame::checked_position_count; +using SocketFrame::make_ipi_virial; +using SocketFrame::transpose_matrix9; +using SocketFrame::validate_ipi_cell; +using SocketFrame::validate_positions; + +const double EPSILON = std::numeric_limits::epsilon(); + +CellValidation validate_with_driver_thresholds(const Matrix9& cell, const Matrix9& inverse) +{ + return validate_ipi_cell(cell, inverse, 1.0e12, 64.0 * EPSILON, 64.0); +} + +void expect_matrix_near(const Matrix9& expected, const Matrix9& actual, double tolerance) +{ + for (std::size_t index = 0; index < expected.size(); ++index) + { + EXPECT_NEAR(expected[index], actual[index], tolerance) << "matrix index " << index; + } +} +} // namespace + +TEST(SocketFrameTest, TransposeKeepsAllNineUniqueEntries) +{ + Matrix9 in = {{1, 2, 3, 4, 5, 6, 7, 8, 9}}; + Matrix9 expected = {{1, 4, 7, 2, 5, 8, 3, 6, 9}}; + EXPECT_EQ(expected, transpose_matrix9(in)); +} + +TEST(SocketFrameTest, VirialUsesPositiveHalfVolumeAndWireTranspose) +{ + Matrix9 stress = {{1, 2, 3, 2, 5, 6, 3, 6, 9}}; + VirialConversion out = make_ipi_virial(stress, 4.0, 1e-12, 1e-12); + Matrix9 expected = {{2, 4, 6, 4, 10, 12, 6, 12, 18}}; + ASSERT_TRUE(out.ok) << out.message; + EXPECT_EQ(expected, out.wire_virial_hartree); +} + +TEST(SocketFrameTest, RightHandedTriclinicCellReturnsKnownInverse) +{ + const Matrix9 cell = {{2.0, 1.0, 0.0, 0.0, 3.0, 1.0, 0.0, 0.0, 4.0}}; + const Matrix9 inverse = {{0.5, -1.0 / 6.0, 1.0 / 24.0, + 0.0, 1.0 / 3.0, -1.0 / 12.0, + 0.0, 0.0, 0.25}}; + + const CellValidation out = validate_with_driver_thresholds(cell, inverse); + + ASSERT_TRUE(out.ok) << out.message; + EXPECT_DOUBLE_EQ(24.0, out.determinant_bohr3); + EXPECT_NEAR(0.0, out.inverse_residual, 16.0 * EPSILON); + expect_matrix_near(inverse, out.computed_inverse_wire_bohr_inv, 16.0 * EPSILON); +} + +TEST(SocketFrameTest, RotatedDiagonalTracksRightSingularVectorsAndInverseOrder) +{ + // Hand-multiplied U diag(5, 2, 0.5) V^T, with rational plane rotations. + const Matrix9 cell = {{4.0, -0.72, -0.96, + 3.0, 0.96, 1.28, + 0.0, -0.4, 0.3}}; + const Matrix9 inverse = {{0.16, 0.12, 0.0, + -0.18, 0.24, -1.6, + -0.24, 0.32, 1.2}}; + + const CellValidation out = validate_with_driver_thresholds(cell, inverse); + + ASSERT_TRUE(out.ok) << out.message; + EXPECT_NEAR(5.0, out.determinant_bohr3, 64.0 * EPSILON); + EXPECT_NEAR(10.0, out.condition_number_2, 256.0 * EPSILON); + expect_matrix_near(inverse, out.computed_inverse_wire_bohr_inv, 64.0 * EPSILON); +} + +TEST(SocketFrameTest, InconsistentReceivedInverseIsRejected) +{ + const Matrix9 cell = {{2.0, 1.0, 0.0, 0.0, 3.0, 1.0, 0.0, 0.0, 4.0}}; + const Matrix9 wrong_inverse = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; + + const CellValidation out = validate_with_driver_thresholds(cell, wrong_inverse); + + EXPECT_FALSE(out.ok); + EXPECT_NE(std::string::npos, out.message.find("inverse")); + EXPECT_GT(out.inverse_residual, 1.0); +} + +TEST(SocketFrameTest, NegativeAndZeroDeterminantsAreRejected) +{ + const Matrix9 identity = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; + const Matrix9 left_handed = {{-1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; + const Matrix9 singular = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}}; + + EXPECT_FALSE(validate_with_driver_thresholds(left_handed, identity).ok); + EXPECT_FALSE(validate_with_driver_thresholds(singular, identity).ok); +} + +TEST(SocketFrameTest, NonrepresentablePositiveCellVolumeIsRejected) +{ + const Matrix9 huge_cell = {{1.0e200, 0.0, 0.0, + 0.0, 1.0e200, 0.0, + 0.0, 0.0, 1.0e200}}; + const Matrix9 tiny_inverse = {{1.0e-200, 0.0, 0.0, + 0.0, 1.0e-200, 0.0, + 0.0, 0.0, 1.0e-200}}; + + const CellValidation out = validate_with_driver_thresholds(huge_cell, tiny_inverse); + + EXPECT_FALSE(out.ok); + EXPECT_NE(std::string::npos, out.message.find("determinant")); +} + +TEST(SocketFrameTest, UnderflowedCellVolumeIsRejectedAsZeroDeterminant) +{ + const Matrix9 tiny_cell = {{1.0e-200, 0.0, 0.0, + 0.0, 1.0e-200, 0.0, + 0.0, 0.0, 1.0e-200}}; + const Matrix9 huge_inverse = {{1.0e200, 0.0, 0.0, + 0.0, 1.0e200, 0.0, + 0.0, 0.0, 1.0e200}}; + + const CellValidation out = validate_with_driver_thresholds(tiny_cell, huge_inverse); + + EXPECT_FALSE(out.ok); + EXPECT_NE(std::string::npos, out.message.find("determinant")); +} + +TEST(SocketFrameTest, ConditionNumberBoundaryIsInclusive) +{ + const Matrix9 below = {{1.0, 0.0, 0.0, 0.0, 1.0e-6, 0.0, 0.0, 0.0, 2.0e-12}}; + const Matrix9 below_inverse = {{1.0, 0.0, 0.0, 0.0, 1.0e6, 0.0, 0.0, 0.0, 5.0e11}}; + const Matrix9 at = {{1.0, 0.0, 0.0, 0.0, 1.0e-6, 0.0, 0.0, 0.0, 1.0e-12}}; + const Matrix9 at_inverse = {{1.0, 0.0, 0.0, 0.0, 1.0e6, 0.0, 0.0, 0.0, 1.0e12}}; + const Matrix9 above = {{1.0, 0.0, 0.0, 0.0, 1.0e-6, 0.0, 0.0, 0.0, 5.0e-13}}; + const Matrix9 above_inverse = {{1.0, 0.0, 0.0, 0.0, 1.0e6, 0.0, 0.0, 0.0, 2.0e12}}; + + EXPECT_TRUE(validate_with_driver_thresholds(below, below_inverse).ok); + const CellValidation boundary = validate_with_driver_thresholds(at, at_inverse); + EXPECT_TRUE(boundary.ok) << boundary.message; + EXPECT_DOUBLE_EQ(1.0e12, boundary.condition_number_2); + EXPECT_FALSE(validate_with_driver_thresholds(above, above_inverse).ok); +} + +TEST(SocketFrameTest, NonfiniteCellOrReceivedInverseIsRejected) +{ + const double nan = std::numeric_limits::quiet_NaN(); + const double infinity = std::numeric_limits::infinity(); + const Matrix9 identity = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; + Matrix9 bad_cell = identity; + Matrix9 bad_inverse = identity; + bad_cell[4] = nan; + EXPECT_FALSE(validate_with_driver_thresholds(bad_cell, identity).ok); + bad_cell = identity; + bad_cell[7] = infinity; + EXPECT_FALSE(validate_with_driver_thresholds(bad_cell, identity).ok); + bad_inverse[1] = nan; + EXPECT_FALSE(validate_with_driver_thresholds(identity, bad_inverse).ok); + bad_inverse = identity; + bad_inverse[8] = -infinity; + EXPECT_FALSE(validate_with_driver_thresholds(identity, bad_inverse).ok); +} + +TEST(SocketFrameTest, PositionCountRequiresMatchingNonnegativeAtomCount) +{ + std::size_t coordinate_count = 77; + std::string message; + + EXPECT_FALSE(checked_position_count(-1, 2, coordinate_count, message)); + EXPECT_EQ(77u, coordinate_count); + EXPECT_NE(std::string::npos, message.find("match")); + + message.clear(); + EXPECT_FALSE(checked_position_count(-1, -1, coordinate_count, message)); + EXPECT_EQ(77u, coordinate_count); + EXPECT_NE(std::string::npos, message.find("negative")); + + message.clear(); + EXPECT_TRUE(checked_position_count(3, 3, coordinate_count, message)) << message; + EXPECT_EQ(9u, coordinate_count); + EXPECT_TRUE(message.empty()); +} + +TEST(SocketFrameTest, PositionCountRejectsMismatchBeforeDerivingAllocationSize) +{ + std::size_t coordinate_count = 123; + std::string message; + + EXPECT_FALSE(checked_position_count(std::numeric_limits::max(), + 1, + coordinate_count, + message)); + EXPECT_EQ(123u, coordinate_count); + EXPECT_NE(std::string::npos, message.find("match")); +} + +TEST(SocketFrameTest, PositionsRequireExactSizeAndFiniteCoordinates) +{ + std::string message; + const std::vector valid = {1.0, -2.0, 3.0}; + EXPECT_TRUE(validate_positions(valid, 3, message)) << message; + + message.clear(); + EXPECT_FALSE(validate_positions(valid, 6, message)); + EXPECT_NE(std::string::npos, message.find("count")); + + std::vector nonfinite = valid; + nonfinite[1] = std::numeric_limits::quiet_NaN(); + message.clear(); + EXPECT_FALSE(validate_positions(nonfinite, 3, message)); + EXPECT_NE(std::string::npos, message.find("finite")); + + nonfinite[1] = std::numeric_limits::infinity(); + message.clear(); + EXPECT_FALSE(validate_positions(nonfinite, 3, message)); + EXPECT_NE(std::string::npos, message.find("finite")); +} + +TEST(SocketFrameTest, SmallStressAsymmetryIsAveragedBeforeConversion) +{ + const Matrix9 stress = {{1.0, 2.1, 3.2, + 1.9, 5.0, 6.3, + 2.8, 5.7, 9.0}}; + const Matrix9 expected = {{1.0, 2.0, 3.0, + 2.0, 5.0, 6.0, + 3.0, 6.0, 9.0}}; + + const VirialConversion out = make_ipi_virial(stress, 2.0, 0.61, 0.0); + + ASSERT_TRUE(out.ok) << out.message; + expect_matrix_near(expected, out.wire_virial_hartree, 4.0 * EPSILON); + EXPECT_NEAR(0.6, out.max_antisymmetric_component, 4.0 * EPSILON); +} + +TEST(SocketFrameTest, ExcessiveStressAsymmetryIsRejected) +{ + const Matrix9 stress = {{1.0, 2.1, 3.2, + 1.9, 5.0, 6.3, + 2.8, 5.7, 9.0}}; + + const VirialConversion out = make_ipi_virial(stress, 2.0, 0.59, 0.0); + + EXPECT_FALSE(out.ok); + EXPECT_NE(std::string::npos, out.message.find("symmetric")); + EXPECT_NEAR(0.6, out.max_antisymmetric_component, 4.0 * EPSILON); +} + +TEST(SocketFrameTest, NonpositiveOrNonfiniteVolumeIsRejected) +{ + const Matrix9 zero_stress = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}; + + EXPECT_FALSE(make_ipi_virial(zero_stress, 0.0, 1.0e-10, 1.0e-8).ok); + EXPECT_FALSE(make_ipi_virial(zero_stress, -1.0, 1.0e-10, 1.0e-8).ok); + EXPECT_FALSE(make_ipi_virial(zero_stress, + std::numeric_limits::infinity(), + 1.0e-10, + 1.0e-8) + .ok); +} + +TEST(SocketFrameTest, NonfiniteStressIsRejected) +{ + Matrix9 stress = {{1.0, 2.0, 3.0, 2.0, 5.0, 6.0, 3.0, 6.0, 9.0}}; + stress[2] = std::numeric_limits::quiet_NaN(); + EXPECT_FALSE(make_ipi_virial(stress, 4.0, 1.0e-10, 1.0e-8).ok); + stress[2] = std::numeric_limits::infinity(); + EXPECT_FALSE(make_ipi_virial(stress, 4.0, 1.0e-10, 1.0e-8).ok); +} From 54ddefdd2edaad1b3477e3338ccebdc217a517c7 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Tue, 4 Aug 2026 00:58:35 +0800 Subject: [PATCH 14/24] fix: enforce socket frame numerical thresholds --- source/source_relax/socket_frame.cpp | 4 +- .../source_relax/test/socket_frame_test.cpp | 67 ++++++++++++++++++- 2 files changed, 67 insertions(+), 4 deletions(-) diff --git a/source/source_relax/socket_frame.cpp b/source/source_relax/socket_frame.cpp index 466499c198e..539cc8bf676 100644 --- a/source/source_relax/socket_frame.cpp +++ b/source/source_relax/socket_frame.cpp @@ -250,9 +250,9 @@ CellValidation validate_ipi_cell(const Matrix9& cell_wire, } result.condition_number_2 = largest_singular / smallest_singular; if (!std::isfinite(result.condition_number_2) - || result.condition_number_2 > max_condition_number) + || result.condition_number_2 >= max_condition_number) { - result.message = "cell condition number exceeds the configured maximum"; + result.message = "cell condition number is not below the configured maximum"; return result; } diff --git a/source/source_relax/test/socket_frame_test.cpp b/source/source_relax/test/socket_frame_test.cpp index 26e25938e58..bba3cef3f6e 100644 --- a/source/source_relax/test/socket_frame_test.cpp +++ b/source/source_relax/test/socket_frame_test.cpp @@ -96,6 +96,30 @@ TEST(SocketFrameTest, InconsistentReceivedInverseIsRejected) EXPECT_GT(out.inverse_residual, 1.0); } +TEST(SocketFrameTest, ReceivedInverseResidualUsesConditionScaledRelativeTolerance) +{ + const Matrix9 cell = {{1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0e-6}}; + Matrix9 accepted_inverse = {{1.0 + 1.0e-8, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0e6}}; + Matrix9 rejected_inverse = accepted_inverse; + rejected_inverse[0] = 1.0 + 2.0e-8; + + const CellValidation accepted + = validate_ipi_cell(cell, accepted_inverse, 1.0e12, 0.0, 64.0); + const CellValidation rejected + = validate_ipi_cell(cell, rejected_inverse, 1.0e12, 0.0, 64.0); + + ASSERT_TRUE(accepted.ok) << accepted.message; + EXPECT_DOUBLE_EQ(1.0e6, accepted.condition_number_2); + EXPECT_NEAR(1.0e-8, accepted.inverse_residual, EPSILON); + EXPECT_FALSE(rejected.ok); + EXPECT_NE(std::string::npos, rejected.message.find("inverse")); + EXPECT_NEAR(2.0e-8, rejected.inverse_residual, EPSILON); +} + TEST(SocketFrameTest, NegativeAndZeroDeterminantsAreRejected) { const Matrix9 identity = {{1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}}; @@ -136,7 +160,7 @@ TEST(SocketFrameTest, UnderflowedCellVolumeIsRejectedAsZeroDeterminant) EXPECT_NE(std::string::npos, out.message.find("determinant")); } -TEST(SocketFrameTest, ConditionNumberBoundaryIsInclusive) +TEST(SocketFrameTest, ConditionNumberMustBeStrictlyBelowMaximum) { const Matrix9 below = {{1.0, 0.0, 0.0, 0.0, 1.0e-6, 0.0, 0.0, 0.0, 2.0e-12}}; const Matrix9 below_inverse = {{1.0, 0.0, 0.0, 0.0, 1.0e6, 0.0, 0.0, 0.0, 5.0e11}}; @@ -147,7 +171,8 @@ TEST(SocketFrameTest, ConditionNumberBoundaryIsInclusive) EXPECT_TRUE(validate_with_driver_thresholds(below, below_inverse).ok); const CellValidation boundary = validate_with_driver_thresholds(at, at_inverse); - EXPECT_TRUE(boundary.ok) << boundary.message; + EXPECT_FALSE(boundary.ok); + EXPECT_NE(std::string::npos, boundary.message.find("condition")); EXPECT_DOUBLE_EQ(1.0e12, boundary.condition_number_2); EXPECT_FALSE(validate_with_driver_thresholds(above, above_inverse).ok); } @@ -255,6 +280,31 @@ TEST(SocketFrameTest, ExcessiveStressAsymmetryIsRejected) EXPECT_NEAR(0.6, out.max_antisymmetric_component, 4.0 * EPSILON); } +TEST(SocketFrameTest, StressAsymmetryUsesAbsolutePlusRelativeTolerance) +{ + const Matrix9 accepted_stress = {{10.0, 2.0 + 4.0e-8, 3.0, + 2.0 - 4.0e-8, 5.0, 6.0, + 3.0, 6.0, 9.0}}; + Matrix9 rejected_stress = accepted_stress; + rejected_stress[1] = 2.0 + 6.0e-8; + rejected_stress[3] = 2.0 - 6.0e-8; + const Matrix9 expected = {{10.0, 2.0, 3.0, + 2.0, 5.0, 6.0, + 3.0, 6.0, 9.0}}; + + const VirialConversion accepted + = make_ipi_virial(accepted_stress, 2.0, 1.0e-10, 1.0e-8); + const VirialConversion rejected + = make_ipi_virial(rejected_stress, 2.0, 1.0e-10, 1.0e-8); + + ASSERT_TRUE(accepted.ok) << accepted.message; + expect_matrix_near(expected, accepted.wire_virial_hartree, 4.0 * EPSILON); + EXPECT_NEAR(8.0e-8, accepted.max_antisymmetric_component, EPSILON); + EXPECT_FALSE(rejected.ok); + EXPECT_NE(std::string::npos, rejected.message.find("symmetric")); + EXPECT_NEAR(1.2e-7, rejected.max_antisymmetric_component, EPSILON); +} + TEST(SocketFrameTest, NonpositiveOrNonfiniteVolumeIsRejected) { const Matrix9 zero_stress = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}}; @@ -268,6 +318,19 @@ TEST(SocketFrameTest, NonpositiveOrNonfiniteVolumeIsRejected) .ok); } +TEST(SocketFrameTest, FiniteStressAndVolumeRejectConvertedVirialOverflow) +{ + const double largest_finite = std::numeric_limits::max(); + const Matrix9 stress = {{largest_finite, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0}}; + + const VirialConversion out = make_ipi_virial(stress, 4.0, 1.0e-10, 1.0e-8); + + EXPECT_FALSE(out.ok); + EXPECT_NE(std::string::npos, out.message.find("representable")); +} + TEST(SocketFrameTest, NonfiniteStressIsRejected) { Matrix9 stress = {{1.0, 2.0, 3.0, 2.0, 5.0, 6.0, 3.0, 6.0, 9.0}}; From 871efbab30002c0949b7f028bd065d31b4e79099 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Wed, 12 Aug 2026 19:46:04 +0800 Subject: [PATCH 15/24] test: cover i-PI extras byte payloads --- source/source_relax/socket_ipi.cpp | 8 ++++++ source/source_relax/socket_ipi.h | 1 + source/source_relax/test/socket_ipi_test.cpp | 29 ++++++++++++++++++++ 3 files changed, 38 insertions(+) diff --git a/source/source_relax/socket_ipi.cpp b/source/source_relax/socket_ipi.cpp index 1c2788dae6c..d0ba19869aa 100644 --- a/source/source_relax/socket_ipi.cpp +++ b/source/source_relax/socket_ipi.cpp @@ -235,6 +235,14 @@ std::string IpiSocket::read_string(std::size_t nbytes) return value; } +void IpiSocket::write_string(const std::string& value) +{ + if (!value.empty()) + { + this->write_exact(value.data(), value.size()); + } +} + void IpiSocket::read_exact(void* data, std::size_t nbytes) { char* cursor = static_cast(data); diff --git a/source/source_relax/socket_ipi.h b/source/source_relax/socket_ipi.h index 563a111f95c..ea8183fd466 100644 --- a/source/source_relax/socket_ipi.h +++ b/source/source_relax/socket_ipi.h @@ -37,6 +37,7 @@ class IpiSocket std::vector read_doubles(std::size_t n); void write_doubles(const std::vector& values); std::string read_string(std::size_t nbytes); + void write_string(const std::string& value); private: int fd_ = -1; diff --git a/source/source_relax/test/socket_ipi_test.cpp b/source/source_relax/test/socket_ipi_test.cpp index cbe92707560..c2b7f415334 100644 --- a/source/source_relax/test/socket_ipi_test.cpp +++ b/source/source_relax/test/socket_ipi_test.cpp @@ -443,3 +443,32 @@ TEST(IpiSocketTest, ReadDoublesRejectsByteCountOverflow) FAIL() << "overflowing double payload size should throw std::overflow_error"; } } + +TEST(IpiSocketTest, WriteStringSendsExactBytesWithoutTerminator) +{ + UnixSocketServer server; + const std::string expected = "{\"scf_converged\":false}"; + std::vector received(expected.size()); + std::exception_ptr thread_error; + std::thread peer([&]() { + try + { + const int fd = server.accept_once(); + recv_all(fd, received.data(), received.size()); + ::close(fd); + } + catch (...) + { + thread_error = std::current_exception(); + } + }); + + IpiSocket socket; + socket.connect(server.address()); + socket.write_string(expected); + socket.close(); + + peer.join(); + rethrow_thread_error(thread_error); + EXPECT_EQ(expected, std::string(received.begin(), received.end())); +} From e16ca567abdd3b660f80cba169e24bd71bca97a5 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Wed, 12 Aug 2026 20:22:12 +0800 Subject: [PATCH 16/24] feat: support independent socket properties --- docs/advanced/input_files/input-main.md | 3 + docs/advanced/interface/ase.md | 6 +- docs/parameters.yaml | 4 + interfaces/ASE_interface/abacuslite/core.py | 253 ++++++++- .../read_input_item_system.cpp | 2 +- source/source_relax/socket_driver.cpp | 494 +++++++++++++--- source/source_relax/test/CMakeLists.txt | 12 + .../source_relax/test/socket_driver_test.cpp | 528 ++++++++++++++++++ 8 files changed, 1204 insertions(+), 98 deletions(-) create mode 100644 source/source_relax/test/socket_driver_test.cpp diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index c33001d4339..36d7d9d8455 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -604,6 +604,7 @@ - host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. - path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. +- Socket mode always computes energy. Force and stress extraction follow `cal_force` and `cal_stress` independently. A disabled property is sent only as fixed-protocol padding and is marked absent in i-PI extras metadata, so a client must not interpret that padding as a physical zero. If SCF does not converge, ABACUS returns the available frame with `scf_converged=false` metadata; the external optimizer or integrator decides whether to accept it. - **Default**: False ### esolver_type @@ -654,6 +655,7 @@ - **Type**: Boolean - **Description**: If set to True, calculate the force at the end of the electronic iteration. +- In `socket_driver` mode, this controls whether forces are calculated and advertised; it is not forced on by the socket protocol. - **Default**: False ### kpar @@ -756,6 +758,7 @@ - **Type**: Boolean - **Description**: If set to True, calculate the stress at the end of the electronic iteration. +- In `socket_driver` mode, this independently controls whether stress/virial is calculated and advertised. - **Default**: False ### diago_proc diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index b9735279d1c..57a3f5841e5 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -105,7 +105,7 @@ Please read the examples in `interfaces/ASE_interface/examples/` for more detail ### Socket I/O with ASE -For socket-driven ASE workflows, use the `AbacusSocketIO` calculator. ASE runs the i-PI socket server, while ABACUS keeps `calculation=scf` and is launched with `socket_driver=1` as the client. The protocol is simple: ASE sends atomic positions and cell data to ABACUS; ABACUS evaluates one SCF step for that structure and returns energy, forces, and virial. See the [ASE socket I/O documentation](https://ase-lib.org/ase/calculators/socketio/socketio.html) and the i-PI reference paper, [Ceriotti et al., Comput. Phys. Commun. 185, 1019-1026 (2014)](https://doi.org/10.1016/j.cpc.2013.10.027), for the protocol background. +For socket-driven ASE workflows, use the `AbacusSocketIO` calculator. ASE runs the i-PI socket server, while ABACUS keeps `calculation=scf` and is launched with `socket_driver=1` as the client. Energy, forces, and stress are independent properties controlled by `cal_force` and `cal_stress`; the fixed i-PI wire layout still contains padding fields, while extras metadata identifies which values were actually computed. See the [ASE socket I/O documentation](https://ase-lib.org/ase/calculators/socketio/socketio.html) and the i-PI reference paper, [Ceriotti et al., Comput. Phys. Commun. 185, 1019-1026 (2014)](https://doi.org/10.1016/j.cpc.2013.10.027), for the protocol background. Build ABACUS as usual before using this interface. PW-only builds work with `basis_type=pw`; LCAO socket calculations require an LCAO-enabled executable. No extra socket library is required. @@ -157,14 +157,14 @@ with abacus as calc: BFGS(atoms).run(fmax=0.05) ``` -`AbacusSocketIO` sets `socket_driver=1` and `cal_force=1` automatically. It also selects the socket endpoint and passes it to ABACUS through `ABACUS_SOCKET_ADDRESS`, so users normally do not set this environment variable by hand when using abacuslite. +`AbacusSocketIO` sets `socket_driver=1` automatically. Set `inp={'cal_force': 1}` and/or `inp={'cal_stress': 1}` when an external optimizer, MD integrator, UnitCellFilter/FrechetCellFilter, or other client needs those properties. Energy is always available. The interface selects the socket endpoint and passes it to ABACUS through `ABACUS_SOCKET_ADDRESS`, so users normally do not set this environment variable by hand when using abacuslite. There are two endpoint styles: - `unixsocket="abacus_si"` uses a local Unix-domain socket. ASE creates and listens on `/tmp/ipi_abacus_si`; abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=/tmp/ipi_abacus_si:UNIX`. The `:UNIX` suffix is part of ABACUS' address syntax and means that `/tmp/ipi_abacus_si` is a filesystem socket path, not a TCP host. This is usually the best choice when ASE and ABACUS run on the same node because it avoids TCP port conflicts. - `port=31415` uses a TCP socket. abacuslite launches ABACUS with `ABACUS_SOCKET_ADDRESS=localhost:31415`, meaning host `localhost` and TCP port `31415`. Use this style when the socket server should listen on a TCP port. If ABACUS is launched manually instead of through `AbacusSocketIO`, set `ABACUS_SOCKET_ADDRESS` yourself to the same `host:port` or `path:UNIX` endpoint. -Calling `atoms.get_potential_energy()` is supported, but the ABACUS client still computes forces because the i-PI `GETFORCE` exchange returns energy, forces, and virial as one response. +Calling `atoms.get_potential_energy()` does not force a force or stress calculation. If a requested property was disabled, ASE raises `PropertyNotImplementedError`; zero-filled i-PI padding is never treated as a physical result. When SCF does not converge, `AbacusSocketIO.last_scf_converged` is set to `False` and the caller decides whether to continue or stop. A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 3cc5af44990..138832d008d 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -48,6 +48,8 @@ parameters: * host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. * path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. + + Socket mode always computes energy. Force and stress extraction follows cal_force and cal_stress independently; disabled properties are sent as protocol padding and marked absent in the i-PI extras metadata, not reported as physical zero values. A non-converged SCF step is returned with scf_converged=false metadata so an external driver can choose its policy. default_value: "False" unit: "" availability: "" @@ -106,6 +108,7 @@ parameters: type: Boolean description: | If set to True, calculate the force at the end of the electronic iteration. + In socket_driver mode, this flag controls whether the returned frame advertises forces; it is not forced on by the socket protocol. default_value: "False" unit: "" availability: "" @@ -222,6 +225,7 @@ parameters: type: Boolean description: | If set to True, calculate the stress at the end of the electronic iteration. + In socket_driver mode, this flag independently controls whether the returned frame advertises stress/virial. default_value: "False" unit: "" availability: "" diff --git a/interfaces/ASE_interface/abacuslite/core.py b/interfaces/ASE_interface/abacuslite/core.py index ca3ead454a7..e1cae67fc82 100644 --- a/interfaces/ASE_interface/abacuslite/core.py +++ b/interfaces/ASE_interface/abacuslite/core.py @@ -30,6 +30,7 @@ @author: Huang Yi-ke ''' +import json import os import re import shutil @@ -594,9 +595,10 @@ class AbacusSocketIO(SocketIOCalculator): A socket calculator owns one ABACUS process with one fixed INPUT/STRU setup. The i-PI protocol can update positions, but electronic-structure parameters such as k-points, spin, basis, pseudopotentials, and species - require a new calculator instance. Energy-only ASE calls are accepted, but - ABACUS still computes forces because i-PI GETFORCE returns energy, forces, - and virial together. + require a new calculator instance. Energy, forces, and stress are + independently controlled by ABACUS INPUT. The fixed-layout i-PI response + uses zero padding for absent fields and an extras metadata record so + padding is never exposed as a computed property. """ def __init__(self, @@ -607,7 +609,19 @@ def __init__(self, timeout=None, log=None, **kwargs): - inp = self._socket_inp(kwargs.pop('inp', {})) + inp = dict(kwargs.pop('inp', {})) + self._property_constraints = {} + for keyword, property_name in (('cal_force', 'forces'), + ('cal_stress', 'stress')): + if keyword in inp: + self._property_constraints[property_name] = self._input_bool( + inp[keyword], keyword) + self.implemented_properties = [ + 'energy', 'free_energy', 'forces', 'stress'] + self._active_properties = None + self._last_socket_metadata = None + self.last_scf_converged = None + inp = self._socket_inp(inp) self.abacus = Abacus( profile=profile, directory=directory, @@ -623,7 +637,7 @@ def __init__(self, launch_client=self._launch_client, ) - def calculate(self, atoms=None, properties=['energy'], system_changes=None): + def calculate(self, atoms=None, properties=None, system_changes=None): from ase.calculators.calculator import ( PropertyNotImplementedError, all_changes, @@ -637,6 +651,9 @@ def calculate(self, atoms=None, properties=['energy'], system_changes=None): if atoms is None: raise ValueError('AbacusSocketIO.calculate requires atoms') + requested = self._normalize_socket_properties(properties) + self._check_requested_properties(requested) + bad = [change for change in system_changes if change not in self.supported_changes] if self.atoms is not None and any(bad): @@ -645,6 +662,22 @@ def calculate(self, atoms=None, properties=['energy'], system_changes=None): 'Please create new socket calculator.' .format(bad if len(bad) > 1 else bad[0])) + desired = set(requested) + desired.discard('free_energy') + desired.add('energy') + for property_name, enabled in self._property_constraints.items(): + if enabled: + desired.add(property_name) + active = set(self._active_properties or ()) + if not active: + active.update(desired) + elif not desired.issubset(active): + active.update(desired) + if self.server is not None: + self._close_socket_session() + self._active_properties = tuple( + name for name in ('energy', 'forces', 'stress') if name in active) + self._check_fixed_cell(atoms) order = self._socket_sort_indices(atoms) socket_atoms = atoms[order] @@ -652,21 +685,66 @@ def calculate(self, atoms=None, properties=['energy'], system_changes=None): if self.server is None: self.server = self.launch_server() - proc = self.launch_client(socket_atoms, properties, + proc = self.launch_client(socket_atoms, list(self._active_properties), port=self._port, unixsocket=self._unixsocket) self.server.proc = proc - results = self.server.calculate(socket_atoms) - results['free_energy'] = results['energy'] - virial = results.pop('virial') - if self.atoms.cell.rank == 3 and any(self.atoms.pbc): - vol = atoms.get_volume() - results['stress'] = -full_3x3_to_voigt_6_stress(virial) / vol - if 'forces' in results: - results['forces'] = self._forces_to_input_order( - results['forces'], order) - self.results.update(results) + raw_results = self.server.calculate(socket_atoms) + if not isinstance(raw_results, dict): + raise ValueError('ABACUS socket server returned a non-mapping result') + results = dict(raw_results) + metadata = self._decode_socket_metadata(results.pop('morebytes', None)) + self._last_socket_metadata = metadata + if metadata is None: + present = set(self._active_properties) + self.last_scf_converged = None + else: + present = set(metadata['present']) + self.last_scf_converged = metadata['scf_converged'] + + if 'energy' not in present or 'energy' not in results: + raise ValueError('ABACUS socket response did not provide energy') + energy = float(results['energy']) + if not np.isfinite(energy): + raise ValueError('ABACUS socket energy is not finite') + free_energy = float(results.get('free_energy', energy)) + if not np.isfinite(free_energy): + raise ValueError('ABACUS socket free energy is not finite') + current = {'energy': energy, 'free_energy': free_energy} + + if 'forces' in present: + if 'forces' not in results: + raise ValueError( + 'ABACUS socket metadata advertises forces, but wire response omitted them') + forces = np.asarray(results['forces'], dtype=np.float64) + expected_shape = (len(socket_atoms), 3) + if forces.shape != expected_shape or not np.all(np.isfinite(forces)): + raise ValueError('ABACUS socket forces have invalid shape or values') + current['forces'] = self._forces_to_input_order(forces, order) + + if 'stress' in present: + virial = results.get('virial') + if virial is None: + raise ValueError( + 'ABACUS socket metadata advertises stress, but wire response omitted virial') + if self.atoms.cell.rank != 3 or not any(self.atoms.pbc): + raise PropertyNotImplementedError( + 'ABACUS socket stress requires a periodic rank-3 cell') + virial = np.asarray(virial, dtype=np.float64) + if virial.shape != (3, 3) or not np.all(np.isfinite(virial)): + raise ValueError('ABACUS socket virial is not a finite 3x3 matrix') + vol = float(atoms.get_volume()) + if not np.isfinite(vol) or vol <= 0.0: + raise ValueError('ABACUS socket stress requires a positive cell volume') + current['stress'] = -full_3x3_to_voigt_6_stress(virial) / vol + + self.results = current + missing = [name for name in requested if name not in current] + if missing: + raise PropertyNotImplementedError( + 'ABACUS socket response did not provide requested {}'.format( + ', '.join(missing))) def _check_fixed_cell(self, atoms): from ase.calculators.calculator import PropertyNotImplementedError @@ -693,11 +771,95 @@ def set(self, **kwargs): ) return super().set(**kwargs) + def _check_requested_properties(self, requested): + from ase.calculators.calculator import PropertyNotImplementedError + + constraints = self._property_constraints + keywords = {'forces': 'cal_force', 'stress': 'cal_stress'} + for property_name, keyword in keywords.items(): + if property_name in requested and constraints.get(property_name) is False: + raise PropertyNotImplementedError( + '{}=0 disables requested {}'.format(keyword, property_name)) + + @staticmethod + def _normalize_socket_properties(properties): + from ase.calculators.calculator import PropertyNotImplementedError + + if properties is None: + names = ['energy'] + elif isinstance(properties, str): + names = [properties] + else: + names = list(properties) + if not names: + names = ['energy'] + allowed = {'energy', 'free_energy', 'forces', 'stress'} + unknown = [name for name in names if name not in allowed] + if unknown: + raise PropertyNotImplementedError( + 'ABACUS socket does not implement {}'.format(', '.join(unknown))) + return tuple(dict.fromkeys(names)) + + def _close_socket_session(self): + server = getattr(self, 'server', None) + if server is not None: + close = getattr(server, 'close', None) + if callable(close): + close() + self.server = None + self.results = {} + + @staticmethod + def _decode_socket_metadata(raw): + if raw is None: + return None + if isinstance(raw, str): + payload = raw.encode('utf-8') + elif isinstance(raw, (bytes, bytearray, memoryview)): + payload = bytes(raw) + else: + payload = np.asarray(raw, dtype=np.uint8).tobytes() + if not payload: + return None + try: + metadata = json.loads(payload.decode('utf-8')) + except (UnicodeDecodeError, ValueError) as error: + raise ValueError('ABACUS socket extras are not valid UTF-8 JSON') from error + if not isinstance(metadata, dict): + raise ValueError('ABACUS socket extras must be a JSON object') + if metadata.get('schema') != 'abacus.socket.properties.v1': + raise ValueError('unsupported ABACUS socket extras schema') + present = metadata.get('present') + if not isinstance(present, list): + raise ValueError('ABACUS socket extras present must be a list') + allowed = {'energy', 'forces', 'stress'} + if any(not isinstance(name, str) or name not in allowed for name in present): + raise ValueError('ABACUS socket extras contain an unknown property') + if 'energy' not in present: + raise ValueError('ABACUS socket extras must include energy') + scf_converged = metadata.get('scf_converged') + if not isinstance(scf_converged, bool): + raise ValueError('ABACUS socket extras scf_converged must be Boolean') + return { + 'present': tuple(dict.fromkeys(present)), + 'scf_converged': scf_converged, + } + def _launch_client(self, atoms, properties=None, port=None, unixsocket=None): from subprocess import Popen if properties is None: - properties = self.abacus.template.implemented_properties + properties = list(self._active_properties or ('energy',)) + properties = set(properties) + properties.discard('free_energy') + properties.add('energy') + # The i-PI response has fixed force/virial fields, but ABACUS must be + # told explicitly which expensive quantities to evaluate. Keep these + # switches synchronized with the session mask before writing INPUT. + self.abacus.parameters['cal_force'] = int('forces' in properties) + self.abacus.parameters['cal_stress'] = int('stress' in properties) + properties = [name for name in ('energy', 'forces', 'stress') + if name in properties] directory = Path(self.abacus.directory) directory.mkdir(exist_ok=True, parents=True) @@ -727,13 +889,29 @@ def _socket_inp(inp): calculation = inp.get('calculation', 'scf') if calculation != 'scf': raise ValueError('ABACUS socket I/O requires calculation="scf"') + for keyword in ('cal_force', 'cal_stress'): + if keyword in inp: + inp[keyword] = int(AbacusSocketIO._input_bool(inp[keyword], keyword)) inp.update({ 'calculation': 'scf', 'socket_driver': 1, - 'cal_force': 1, }) return inp + @staticmethod + def _input_bool(value, name): + if isinstance(value, bool): + return value + if isinstance(value, int) and value in (0, 1): + return bool(value) + if isinstance(value, str): + normalized = value.strip().lower() + if normalized in ('true', '1'): + return True + if normalized in ('false', '0'): + return False + raise ValueError('{} must be one of true, false, 1, or 0'.format(name)) + @staticmethod def _socket_sort_indices(atoms): return species_group_indices(atoms.get_chemical_symbols()) @@ -782,6 +960,45 @@ def test_socketio_rejects_cell_changes(self): with self.assertRaisesRegex(PropertyNotImplementedError, 'fixed-cell'): calc._check_fixed_cell(changed) + def test_socketio_input_keeps_independent_property_switches(self): + self.assertEqual( + AbacusSocketIO._socket_inp({'cal_force': 0, 'cal_stress': 1}), + {'calculation': 'scf', 'socket_driver': 1, + 'cal_force': 0, 'cal_stress': 1}) + + def test_socketio_boolean_parser_rejects_ambiguous_values(self): + for value in ('yes', 'no', 'on', 'off', ''): + with self.assertRaisesRegex(ValueError, 'cal_force'): + AbacusSocketIO._input_bool(value, 'cal_force') + + def test_socketio_metadata_rejects_unknown_property(self): + metadata = json.dumps({ + 'schema': 'abacus.socket.properties.v1', + 'present': ['energy', 'charges'], + 'scf_converged': True, + }).encode('utf-8') + with self.assertRaisesRegex(ValueError, 'unknown property'): + AbacusSocketIO._decode_socket_metadata( + np.frombuffer(metadata, dtype=np.int8)) + + def test_socketio_metadata_marks_padding_absent(self): + metadata = json.dumps({ + 'schema': 'abacus.socket.properties.v1', + 'present': ['energy'], + 'scf_converged': False, + }).encode('utf-8') + decoded = AbacusSocketIO._decode_socket_metadata( + np.frombuffer(metadata, dtype=np.int8)) + self.assertEqual(decoded['present'], ('energy',)) + self.assertFalse(decoded['scf_converged']) + + def test_socketio_requested_disabled_property_is_rejected(self): + calc = object.__new__(AbacusSocketIO) + calc._property_constraints = {'forces': False, 'stress': False} + from ase.calculators.calculator import PropertyNotImplementedError + with self.assertRaises(PropertyNotImplementedError): + calc._check_requested_properties(('forces',)) + def test_parse_version_allows_launcher_noise(self): stdout = 'launcher info\nABACUS version v3.11.0-beta6\n' self.assertEqual(AbacusProfile.parse_version(stdout), 'v3.11.0-beta6') diff --git a/source/source_io/module_parameter/read_input_item_system.cpp b/source/source_io/module_parameter/read_input_item_system.cpp index 3fd78ed0d97..e7fe585f275 100644 --- a/source/source_io/module_parameter/read_input_item_system.cpp +++ b/source/source_io/module_parameter/read_input_item_system.cpp @@ -293,7 +293,7 @@ When using the ASE AbacusSocketIO interface, this environment variable is set au item.reset_value = [](const Input_Item& item, Parameter& para) { std::vector use_force = {"cell-relax", "relax", "md"}; std::vector not_use_force = {"get_wf", "get_pchg", "get_s"}; - if (para.input.socket_driver || std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) + if (std::find(use_force.begin(), use_force.end(), para.input.calculation) != use_force.end()) { if (!para.input.cal_force) { diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index b97e4c118d2..2ca114dbd1f 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -1,6 +1,7 @@ #include "socket_driver.h" #include "source_relax/socket_ipi.h" +#include "source_relax/socket_frame.h" #include "source_base/global_function.h" #include "source_base/mathzone.h" #include "source_base/parallel_common.h" @@ -12,8 +13,13 @@ #include #include +#include +#include #include #include +#include +#include +#include #include #include @@ -21,6 +27,93 @@ namespace { constexpr double RY_TO_HARTREE = 0.5; constexpr int IPI_RANK_ROOT = 0; +constexpr double MAX_CELL_CONDITION = 1.0e12; +constexpr double INVERSE_ABSOLUTE_TOLERANCE + = 64.0 * std::numeric_limits::epsilon(); +constexpr double INVERSE_RELATIVE_TOLERANCE = 64.0; +constexpr double STRESS_ABSOLUTE_TOLERANCE = 1.0e-10; +constexpr double STRESS_RELATIVE_TOLERANCE = 1.0e-8; +constexpr std::int32_t MAX_INIT_BYTES = INT32_C(1048576); + +enum class DriverState +{ + NeedInit, + Ready, + HasData +}; + +struct ComputedFrame +{ + bool valid = false; + bool forces_present = false; + bool stress_present = false; + bool scf_converged = true; + double energy_hartree = 0.0; + std::vector forces_hartree_per_bohr; + SocketFrame::Matrix9 virial_wire_hartree = {{0.0}}; +}; + +bool all_ranks_converged(const bool local_converged) +{ + int converged = local_converged ? 1 : 0; +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, &converged, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); +#endif + return converged != 0; +} + +void throw_if_any_rank_failed(int local_failed, std::string local_message) +{ + int any_failed = local_failed; +#ifdef __MPI + MPI_Allreduce(MPI_IN_PLACE, &any_failed, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); +#endif + if (any_failed != 0) + { + if (local_message.empty()) + { + local_message = "socket frame validation failed on another MPI rank"; + } + throw std::runtime_error(local_message); + } +} + +[[noreturn]] void fail_during_collective_stage(const char* stage, + const std::string& message) +{ +#ifdef __MPI + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, + "ABACUS_SOCKET_MPI_FATAL stage=%s rank=%d message=%s\n", + stage, + rank, + message.c_str()); + std::fflush(stderr); + MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE); + std::abort(); +#else + (void)stage; + throw std::runtime_error(message); +#endif +} + +std::string properties_extra(const ComputedFrame& frame) +{ + std::ostringstream extra; + extra << "{\"schema\":\"abacus.socket.properties.v1\",\"present\":[\"energy\""; + if (frame.forces_present) + { + extra << ",\"forces\""; + } + if (frame.stress_present) + { + extra << ",\"stress\""; + } + extra << "],\"scf_converged\":" + << (frame.scf_converged ? "true" : "false") << "}"; + return extra.str(); +} bool is_root() { @@ -54,6 +147,15 @@ void bcast_socket_int(int& value) #endif } +void bcast_socket_int32(std::int32_t& value) +{ +#ifdef __MPI + MPI_Bcast(&value, 1, MPI_INT32_T, IPI_RANK_ROOT, MPI_COMM_WORLD); +#else + (void)value; +#endif +} + void bcast_socket_chars(char* value, const int size) { #ifdef __MPI @@ -182,6 +284,16 @@ double max_abs_delta(const std::vector& a, const std::vector& b) return out; } +double unchanged_cell_tolerance(const SocketFrame::Matrix9& cell) +{ + double maximum = 0.0; + for (std::size_t index = 0; index < cell.size(); ++index) + { + maximum = std::max(maximum, std::fabs(cell[index])); + } + return 32.0 * std::numeric_limits::epsilon() * std::max(1.0, maximum); +} + void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& positions_bohr) { if (positions_bohr.size() != static_cast(3 * ucell.nat)) @@ -237,16 +349,47 @@ void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& pos std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force) { + if (force.nr < 0 || force.nc != 3) + { + throw std::runtime_error("force matrix must have nat rows and three columns"); + } std::vector out(static_cast(force.nr * force.nc)); for (int iat = 0; iat < force.nr; ++iat) { for (int idir = 0; idir < force.nc; ++idir) { - out[static_cast(3 * iat + idir)] = force(iat, idir) * RY_TO_HARTREE; + const double value = force(iat, idir); + if (!std::isfinite(value)) + { + throw std::runtime_error("force entries must be finite"); + } + out[static_cast(3 * iat + idir)] = value * RY_TO_HARTREE; } } return out; } + +SocketFrame::Matrix9 matrix9_from_stress(const ModuleBase::matrix& stress) +{ + if (stress.nr != 3 || stress.nc != 3) + { + throw std::runtime_error("stress matrix must have three rows and three columns"); + } + SocketFrame::Matrix9 values; + for (int row = 0; row < 3; ++row) + { + for (int column = 0; column < 3; ++column) + { + values[3 * row + column] = stress(row, column); + } + } + return values; +} + +std::vector vector_from_matrix9(const SocketFrame::Matrix9& values) +{ + return std::vector(values.begin(), values.end()); +} } // namespace void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, @@ -261,11 +404,6 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { ModuleBase::WARNING_QUIT("ABACUS socket", "socket driver requires a valid ESolver."); } - if (!inp.cal_force) - { - ModuleBase::WARNING_QUIT("ABACUS socket", "socket_driver requires cal_force=1 for i-PI GETFORCE."); - } - IpiSocket socket; try @@ -288,13 +426,10 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } quit_if_root_io_failed(io_failed, io_message); - bool isinit = false; - bool hasdata = false; + DriverState state = DriverState::NeedInit; int istep = 0; const int nat_return = ucell.nat; - double energy_hartree = 0.0; - std::vector forces_hartree_bohr(static_cast(3 * ucell.nat), 0.0); - std::vector virial_hartree(9, 0.0); + ComputedFrame published; const std::vector reference_cell = ipi_cell_bohr_from_unitcell(ucell); bool checked_initial_positions = false; @@ -312,7 +447,15 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } catch (const IpiSocketClosed&) { - header.clear(); + if (state == DriverState::HasData) + { + io_failed = 1; + io_message = "i-PI peer closed while a computed frame was pending"; + } + else + { + header.clear(); + } } catch (const std::exception& exc) { @@ -339,11 +482,11 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { try { - if (hasdata) + if (state == DriverState::HasData) { socket.write_header("HAVEDATA"); } - else if (isinit) + else if (state == DriverState::Ready) { socket.write_header("READY"); } @@ -362,41 +505,54 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } else if (header == "INIT") { - int rid = 0; - int nbytes = 0; + std::int32_t rid = 0; + std::int32_t nbytes = 0; std::string params; io_failed = 0; io_message.clear(); if (is_root()) { - try + if (state != DriverState::NeedInit) + { + io_failed = 1; + io_message = "INIT requires NEEDINIT state"; + } + else { - rid = socket.read_int32(); - nbytes = socket.read_int32(); - if (nbytes < 0) + try { - io_failed = 1; - io_message = "negative INIT payload length from i-PI socket"; + rid = socket.read_int32(); + nbytes = socket.read_int32(); + if (nbytes < 0) + { + io_failed = 1; + io_message = "negative INIT payload length from i-PI socket"; + } + else if (nbytes > MAX_INIT_BYTES) + { + io_failed = 1; + io_message = "INIT payload exceeds the 1 MiB socket limit"; + } + else if (nbytes > 0) + { + params = socket.read_string(static_cast(nbytes)); + } } - else if (nbytes > 0) + catch (const std::exception& exc) { - params = socket.read_string(static_cast(nbytes)); + io_failed = 1; + io_message = exc.what(); } } - catch (const std::exception& exc) - { - io_failed = 1; - io_message = exc.what(); - } } quit_if_root_io_failed(io_failed, io_message); - bcast_socket_int(rid); - bcast_socket_int(nbytes); + bcast_socket_int32(rid); + bcast_socket_int32(nbytes); if (nbytes > 0 && is_root()) { ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; } - isinit = true; + state = DriverState::Ready; if (is_root()) { ofs_running << " ABACUS socket INIT replica " << rid << std::endl; @@ -404,51 +560,87 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } else if (header == "POSDATA") { - std::vector cell(9, 0.0); - std::vector inv_cell(9, 0.0); - int nat_socket = 0; + SocketFrame::Matrix9 cell = {{0.0}}; + SocketFrame::Matrix9 inv_cell = {{0.0}}; + std::int32_t nat_socket = 0; std::vector positions; io_failed = 0; io_message.clear(); if (is_root()) { - try + if (state != DriverState::Ready) + { + io_failed = 1; + io_message = "POSDATA requires READY state"; + } + else { - cell = socket.read_doubles(9); - inv_cell = socket.read_doubles(9); - nat_socket = socket.read_int32(); - if (nat_socket < 0) + try { - io_failed = 1; - io_message = "negative POSDATA atom count from i-PI socket"; + const std::vector cell_values = socket.read_doubles(9); + const std::vector inverse_values = socket.read_doubles(9); + std::copy(cell_values.begin(), cell_values.end(), cell.begin()); + std::copy(inverse_values.begin(), inverse_values.end(), inv_cell.begin()); + nat_socket = socket.read_int32(); + SocketFrame::CellValidation validation + = SocketFrame::validate_ipi_cell(cell, + inv_cell, + MAX_CELL_CONDITION, + INVERSE_ABSOLUTE_TOLERANCE, + INVERSE_RELATIVE_TOLERANCE); + if (!validation.ok) + { + io_failed = 1; + io_message = "invalid POSDATA cell: " + validation.message; + } + std::size_t coordinate_count = 0; + if (io_failed == 0 + && !SocketFrame::checked_position_count(nat_socket, + ucell.nat, + coordinate_count, + io_message)) + { + io_failed = 1; + } + if (io_failed == 0) + { + positions = socket.read_doubles(coordinate_count); + if (!SocketFrame::validate_positions(positions, + coordinate_count, + io_message)) + { + io_failed = 1; + } + } } - else + catch (const std::exception& exc) { - positions = socket.read_doubles(static_cast(3 * nat_socket)); + io_failed = 1; + io_message = exc.what(); } } - catch (const std::exception& exc) - { - io_failed = 1; - io_message = exc.what(); - } } quit_if_root_io_failed(io_failed, io_message); - bcast_double_vector(cell); - bcast_double_vector(inv_cell); - bcast_socket_int(nat_socket); + bcast_socket_int32(nat_socket); + std::vector cell_values(cell.begin(), cell.end()); + std::vector inverse_values(inv_cell.begin(), inv_cell.end()); + bcast_double_vector(cell_values); + bcast_double_vector(inverse_values); if (!is_root()) { - positions.assign(static_cast(3 * nat_socket), 0.0); + cell = {{0.0}}; + inv_cell = {{0.0}}; + std::copy(cell_values.begin(), cell_values.end(), cell.begin()); + std::copy(inverse_values.begin(), inverse_values.end(), inv_cell.begin()); + if (nat_socket >= 0) + { + positions.assign(static_cast(3 * nat_socket), 0.0); + } } bcast_double_vector(positions); - if (nat_socket != ucell.nat) - { - ModuleBase::WARNING_QUIT("ABACUS socket", "POSDATA atom count does not match STRU."); - } - const double max_cell_delta_bohr = max_abs_delta(cell, reference_cell); - if (max_cell_delta_bohr > 1.0e-6) + const double max_cell_delta_bohr = max_abs_delta(std::vector(cell.begin(), cell.end()), reference_cell); + if (max_cell_delta_bohr > unchanged_cell_tolerance(cell)) { ModuleBase::WARNING_QUIT("ABACUS socket", "variable-cell socket updates are not supported yet."); } @@ -464,25 +656,152 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } } - set_positions_from_ipi_bohr(ucell, positions); - p_esolver->runner(ucell, istep); - const double energy_ry = p_esolver->cal_energy(); - energy_hartree = energy_ry * RY_TO_HARTREE; + try + { + set_positions_from_ipi_bohr(ucell, positions); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("set_positions", exc.what()); + } + catch (...) + { + fail_during_collective_stage("set_positions", + "unknown socket position update failure"); + } + try + { + p_esolver->runner(ucell, istep); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("runner", exc.what()); + } + catch (...) + { + fail_during_collective_stage("runner", + "unknown socket runner failure"); + } + ComputedFrame computed; + computed.scf_converged = all_ranks_converged(p_esolver->conv_esolver); + if (!computed.scf_converged && is_root()) + { + ModuleBase::WARNING( + "ABACUS socket", + "SCF did not converge; returning the available frame and marking it in i-PI extras."); + } + double energy_ry = 0.0; + try + { + energy_ry = p_esolver->cal_energy(); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("cal_energy", exc.what()); + } + catch (...) + { + fail_during_collective_stage("cal_energy", + "unknown socket energy failure"); + } + int local_failed = std::isfinite(energy_ry) ? 0 : 1; + throw_if_any_rank_failed(local_failed, + local_failed == 0 ? "" : "socket energy is not finite"); + if (!std::isfinite(energy_ry)) + { + ModuleBase::WARNING_QUIT("ABACUS socket", "socket energy is not finite."); + } + computed.energy_hartree = energy_ry * RY_TO_HARTREE; if (is_root()) { ofs_running << " ABACUS socket return energy " << energy_ry << " Ry, " << energy_ry * ModuleBase::Ry_to_eV << " eV, " - << energy_hartree << " Ha" << std::endl; + << computed.energy_hartree << " Ha" << std::endl; } ModuleBase::matrix force; if (inp.cal_force) { - p_esolver->cal_force(ucell, force); - forces_hartree_bohr = flatten_forces_hartree_per_bohr(force); + try + { + p_esolver->cal_force(ucell, force); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("cal_force", exc.what()); + } + catch (...) + { + fail_during_collective_stage("cal_force", + "unknown socket force failure"); + } + local_failed = 0; + std::string local_message; + try + { + computed.forces_hartree_per_bohr = flatten_forces_hartree_per_bohr(force); + } + catch (const std::exception& exc) + { + local_failed = 1; + local_message = exc.what(); + } + catch (...) + { + local_failed = 1; + local_message = "unknown socket force validation failure"; + } + throw_if_any_rank_failed(local_failed, local_message); + computed.forces_present = true; + } + if (inp.cal_stress) + { + ModuleBase::matrix stress; + try + { + p_esolver->cal_stress(ucell, stress); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("cal_stress", exc.what()); + } + catch (...) + { + fail_during_collective_stage("cal_stress", + "unknown socket stress failure"); + } + local_failed = 0; + std::string local_message; + try + { + const SocketFrame::VirialConversion virial + = SocketFrame::make_ipi_virial(matrix9_from_stress(stress), + ucell.omega, + STRESS_ABSOLUTE_TOLERANCE, + STRESS_RELATIVE_TOLERANCE); + if (!virial.ok) + { + throw std::runtime_error(virial.message); + } + computed.virial_wire_hartree = virial.wire_virial_hartree; + } + catch (const std::exception& exc) + { + local_failed = 1; + local_message = exc.what(); + } + catch (...) + { + local_failed = 1; + local_message = "unknown socket stress validation failure"; + } + throw_if_any_rank_failed(local_failed, local_message); + computed.stress_present = true; } + computed.valid = true; + published = computed; ++istep; - hasdata = true; + state = DriverState::HasData; } else if (header == "GETFORCE") { @@ -492,12 +811,26 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { try { + if (state != DriverState::HasData || !published.valid) + { + throw std::runtime_error("GETFORCE requires HAVEDATA state and a valid frame"); + } socket.write_header("FORCEREADY"); - socket.write_double(energy_hartree); - socket.write_int32(nat_return); - socket.write_doubles(forces_hartree_bohr); - socket.write_doubles(virial_hartree); - socket.write_int32(0); + socket.write_double(published.energy_hartree); + socket.write_int32(static_cast(nat_return)); + const std::vector forces + = published.forces_present + ? published.forces_hartree_per_bohr + : std::vector(static_cast(3 * nat_return), 0.0); + socket.write_doubles(forces); + socket.write_doubles(vector_from_matrix9(published.virial_wire_hartree)); + const std::string extra = properties_extra(published); + if (extra.size() > static_cast(std::numeric_limits::max())) + { + throw std::overflow_error("i-PI extras payload is larger than int32"); + } + socket.write_int32(static_cast(extra.size())); + socket.write_string(extra); } catch (const std::exception& exc) { @@ -506,17 +839,26 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } } quit_if_root_io_failed(io_failed, io_message); - isinit = false; - hasdata = false; + published = ComputedFrame(); + state = DriverState::Ready; } - else + else if (header == "EXIT") { if (is_root()) { - ofs_running << " ABACUS socket driver exiting on header " << header << std::endl; + ofs_running << " ABACUS socket driver received i-PI EXIT" << std::endl; } break; } + else + { + if (is_root()) + { + io_failed = 1; + io_message = "unknown i-PI header: " + header; + } + quit_if_root_io_failed(io_failed, io_message); + } } } catch (const std::exception& exc) diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index 3b50983345a..6fdac704870 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -17,6 +17,18 @@ AddTest( SOURCES socket_frame_test.cpp ../socket_frame.cpp ) +AddTest( + TARGET MODULE_RELAX_socket_driver_test + LIBS base device + SOURCES socket_driver_test.cpp + ../socket_driver.cpp + ../socket_frame.cpp + ../socket_ipi.cpp + ../../source_cell/update_cell.cpp + ../../source_cell/bcast_cell.cpp +) +set_tests_properties(MODULE_RELAX_socket_driver_test PROPERTIES TIMEOUT 15) + AddTest( TARGET MODULE_RELAX_relax_new_line_search LIBS parameter diff --git a/source/source_relax/test/socket_driver_test.cpp b/source/source_relax/test/socket_driver_test.cpp new file mode 100644 index 00000000000..7f8c4255c2a --- /dev/null +++ b/source/source_relax/test/socket_driver_test.cpp @@ -0,0 +1,528 @@ +#include "source_relax/socket_driver.h" + +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include "source_cell/unitcell.h" +#include "source_esolver/esolver.h" +#include "source_io/module_parameter/input_parameter.h" +#include "for_test.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace unitcell +{ +void periodic_boundary_adjustment(Atom* atoms, const ModuleBase::Matrix3& latvec, const int ntype) +{ + for (int it = 0; it < ntype; ++it) + { + for (int ia = 0; ia < atoms[it].na; ++ia) + { + atoms[it].tau[ia] = atoms[it].taud[ia] * latvec; + } + } +} +} // namespace unitcell + +namespace +{ +constexpr std::size_t IPI_HEADER_LEN = 12; + +std::string errno_message(const std::string& prefix) +{ + return prefix + ": " + std::strerror(errno); +} + +void send_all(const int fd, const void* data, const std::size_t nbytes) +{ + const char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { +#ifdef MSG_NOSIGNAL + const int flags = MSG_NOSIGNAL; +#else + const int flags = 0; +#endif + const ssize_t sent = ::send(fd, cursor + done, nbytes - done, flags); + if (sent < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("send failed")); + } + if (sent == 0) + { + throw std::runtime_error("send returned zero"); + } + done += static_cast(sent); + } +} + +template +void send_value(const int fd, const T& value) +{ + send_all(fd, &value, sizeof(value)); +} + +void send_header(const int fd, const std::string& header) +{ + std::string padded = header; + padded.resize(IPI_HEADER_LEN, ' '); + send_all(fd, padded.data(), padded.size()); +} + +bool try_send_status(const int fd) +{ + try + { + send_header(fd, "STATUS"); + return true; + } + catch (const std::runtime_error&) + { + if (errno == EPIPE || errno == ECONNRESET) + { + return false; + } + throw; + } +} + +std::string read_header_or_close(const int fd) +{ + char header[IPI_HEADER_LEN]; + std::size_t done = 0; + while (done < sizeof(header)) + { + const ssize_t received = ::recv(fd, header + done, sizeof(header) - done, 0); + if (received == 0 || (received < 0 && errno == ECONNRESET)) + { + if (done == 0) + { + return ""; + } + throw std::runtime_error("socket closed during response header"); + } + if (received < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("receive failed")); + } + done += static_cast(received); + } + + std::string value(header, sizeof(header)); + while (!value.empty() && value.back() == ' ') + { + value.pop_back(); + } + return value; +} + +class UnixSocketServer +{ + public: + UnixSocketServer() + { + char dir_template[] = "/tmp/abacus_socket_driver_test_XXXXXX"; + char* made_dir = ::mkdtemp(dir_template); + if (made_dir == nullptr) + { + throw std::runtime_error(errno_message("mkdtemp failed")); + } + dir_ = made_dir; + path_ = dir_ + "/ipi.sock"; + + listen_fd_ = ::socket(AF_UNIX, SOCK_STREAM, 0); + if (listen_fd_ < 0) + { + throw std::runtime_error(errno_message("socket failed")); + } + + sockaddr_un address; + std::memset(&address, 0, sizeof(address)); + address.sun_family = AF_UNIX; + std::strncpy(address.sun_path, path_.c_str(), sizeof(address.sun_path) - 1); + if (::bind(listen_fd_, reinterpret_cast(&address), sizeof(address)) != 0) + { + throw std::runtime_error(errno_message("bind failed")); + } + if (::listen(listen_fd_, 1) != 0) + { + throw std::runtime_error(errno_message("listen failed")); + } + } + + ~UnixSocketServer() + { + if (listen_fd_ >= 0) + { + ::close(listen_fd_); + } + if (!path_.empty()) + { + ::unlink(path_.c_str()); + } + if (!dir_.empty()) + { + ::rmdir(dir_.c_str()); + } + } + + UnixSocketServer(const UnixSocketServer&) = delete; + UnixSocketServer& operator=(const UnixSocketServer&) = delete; + + std::string address() const + { + return path_ + ":UNIX"; + } + + int accept_once() const + { + const int fd = ::accept(listen_fd_, nullptr, nullptr); + if (fd < 0) + { + throw std::runtime_error(errno_message("accept failed")); + } + return fd; + } + + private: + int listen_fd_ = -1; + std::string dir_; + std::string path_; +}; + +class FakeESolver : public ModuleESolver::ESolver +{ + public: + explicit FakeESolver(const bool converged) : converged_(converged) + { + } + + void before_all_runners(BaseCell&, const Input_para&) override + { + } + + void runner(BaseCell&, const int) override + { + this->conv_esolver = converged_; + } + + void after_all_runners(BaseCell&) override + { + } + + double cal_energy() override + { + return 4.0; + } + + void cal_force(BaseCell& cell, ModuleBase::matrix& force) override + { + force.create(cell.nat(), 3); + } + + void cal_stress(BaseCell&, ModuleBase::matrix& stress) override + { + stress.create(3, 3); + } + + private: + bool converged_; +}; + +struct DriverResult +{ + int exit_code = -1; + std::string response_header; + std::string diagnostic; +}; + +struct ForceResponse +{ + std::string header; + double energy_hartree = 0.0; + std::int32_t nat = 0; + std::vector forces_hartree_per_bohr; + std::vector virial_wire_hartree; + std::string extra; +}; + +void initialize_one_atom_cell(UnitCell& ucell) +{ + ucell.lat0 = 1.0; + ucell.latvec.Identity(); + ucell.ntype = 1; + ucell.nat = 1; + ucell.atoms[0].na = 1; + ucell.atoms[0].tau.resize(1); + ucell.atoms[0].taud.resize(1); + ucell.atoms[0].dis.resize(1); +} + +void send_fixed_cell_frame(const int fd) +{ + const std::int32_t replica = 0; + const std::int32_t parameter_bytes = 0; + send_header(fd, "INIT"); + send_value(fd, replica); + send_value(fd, parameter_bytes); + + const double identity[9] = {1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0}; + const std::int32_t nat = 1; + const double position[3] = {0.0, 0.0, 0.0}; + send_header(fd, "POSDATA"); + send_all(fd, identity, sizeof(identity)); + send_all(fd, identity, sizeof(identity)); + send_value(fd, nat); + send_all(fd, position, sizeof(position)); +} + +void read_all(const int fd, void* data, const std::size_t nbytes) +{ + char* cursor = static_cast(data); + std::size_t done = 0; + while (done < nbytes) + { + const ssize_t received = ::recv(fd, cursor + done, nbytes - done, 0); + if (received <= 0) + { + throw std::runtime_error("socket closed while reading response"); + } + done += static_cast(received); + } +} + +template +T read_value(const int fd) +{ + T value; + read_all(fd, &value, sizeof(value)); + return value; +} + +std::vector read_doubles(const int fd, const std::size_t count) +{ + std::vector values(count); + if (!values.empty()) + { + read_all(fd, values.data(), values.size() * sizeof(double)); + } + return values; +} + +ForceResponse read_force_response(const int fd) +{ + ForceResponse response; + response.header = read_header_or_close(fd); + if (response.header.empty()) + { + return response; + } + response.energy_hartree = read_value(fd); + response.nat = read_value(fd); + response.forces_hartree_per_bohr + = read_doubles(fd, static_cast(3 * response.nat)); + response.virial_wire_hartree = read_doubles(fd, 9); + const std::int32_t extra_bytes = read_value(fd); + if (extra_bytes < 0) + { + throw std::runtime_error("negative extras length"); + } + response.extra.resize(static_cast(extra_bytes)); + if (!response.extra.empty()) + { + read_all(fd, &response.extra[0], response.extra.size()); + } + return response; +} + +std::string read_pipe(const int fd) +{ + std::string output; + char buffer[512]; + while (true) + { + const ssize_t nread = ::read(fd, buffer, sizeof(buffer)); + if (nread == 0) + { + break; + } + if (nread < 0) + { + if (errno == EINTR) + { + continue; + } + throw std::runtime_error(errno_message("pipe read failed")); + } + output.append(buffer, static_cast(nread)); + } + return output; +} + +DriverResult run_driver_frame(const bool converged, + const bool cal_force, + const bool cal_stress, + const std::function& peer_action) +{ + UnixSocketServer server; + int output_pipe[2]; + if (::pipe(output_pipe) != 0) + { + throw std::runtime_error(errno_message("pipe failed")); + } + + const pid_t child = ::fork(); + if (child < 0) + { + ::close(output_pipe[0]); + ::close(output_pipe[1]); + throw std::runtime_error(errno_message("fork failed")); + } + if (child == 0) + { + ::close(output_pipe[0]); + ::dup2(output_pipe[1], STDOUT_FILENO); + ::dup2(output_pipe[1], STDERR_FILENO); + ::close(output_pipe[1]); + ::setenv("ABACUS_SOCKET_ADDRESS", server.address().c_str(), 1); + + UnitCell ucell; + initialize_one_atom_cell(ucell); + Input_para input; + input.cal_force = cal_force; + input.cal_stress = cal_stress; + FakeESolver solver(converged); + std::ofstream running("/dev/null"); + Socket_Driver driver; + driver.socket_driver(&solver, ucell, input, running); + std::cout.flush(); + std::cerr.flush(); + ::_exit(0); + } + + ::close(output_pipe[1]); + DriverResult result; + std::exception_ptr peer_error; + int peer_fd = -1; + try + { + peer_fd = server.accept_once(); + timeval timeout; + timeout.tv_sec = 5; + timeout.tv_usec = 0; + if (::setsockopt(peer_fd, SOL_SOCKET, SO_RCVTIMEO, &timeout, sizeof(timeout)) != 0) + { + throw std::runtime_error(errno_message("setsockopt failed")); + } + send_fixed_cell_frame(peer_fd); + peer_action(peer_fd); + } + catch (...) + { + peer_error = std::current_exception(); + } + if (peer_fd >= 0) + { + ::close(peer_fd); + } + + result.diagnostic = read_pipe(output_pipe[0]); + ::close(output_pipe[0]); + int status = 0; + while (::waitpid(child, &status, 0) < 0) + { + if (errno != EINTR) + { + throw std::runtime_error(errno_message("waitpid failed")); + } + } + if (WIFEXITED(status)) + { + result.exit_code = WEXITSTATUS(status); + } + + if (peer_error) + { + std::rethrow_exception(peer_error); + } + return result; +} +} // namespace + +TEST(SocketDriverTest, NonconvergedFrameIsPublishedWithMetadata) +{ + ForceResponse response; + const DriverResult result = run_driver_frame( + false, true, false, + [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); + + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(response.extra, testing::HasSubstr("\"scf_converged\":false")); +} + +TEST(SocketDriverTest, EnergyOnlyFrameMarksForceAndStressAbsent) +{ + ForceResponse response; + const DriverResult result = run_driver_frame( + true, false, false, + [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); + + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\"]")); + EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"forces\""))); + EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"stress\""))); + EXPECT_THAT(response.forces_hartree_per_bohr, + testing::ElementsAre(0.0, 0.0, 0.0)); + EXPECT_THAT(response.virial_wire_hartree, + testing::ElementsAre(0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0)); +} + +TEST(SocketDriverTest, EnergyAndStressFrameDoesNotAdvertiseForce) +{ + ForceResponse response; + const DriverResult result = run_driver_frame( + true, false, true, + [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); + + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"stress\"]")); + EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"forces\""))); + EXPECT_NE(0.0, response.virial_wire_hartree[0]); +} From 1e10c5970edf5e462a818b169d56b53853e877e5 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Wed, 12 Aug 2026 20:39:50 +0800 Subject: [PATCH 17/24] test: fix socket driver integration linkage --- source/source_relax/socket_driver.cpp | 6 +++--- .../source_relax/test/socket_driver_test.cpp | 18 ++++-------------- 2 files changed, 7 insertions(+), 17 deletions(-) diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index 2ca114dbd1f..a5bf94cde21 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -347,9 +347,9 @@ void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& pos ucell.cell_parameter_updated = false; } -std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force) +std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force, const int nat) { - if (force.nr < 0 || force.nc != 3) + if (nat < 0 || force.nr != nat || force.nc != 3) { throw std::runtime_error("force matrix must have nat rows and three columns"); } @@ -739,7 +739,7 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, std::string local_message; try { - computed.forces_hartree_per_bohr = flatten_forces_hartree_per_bohr(force); + computed.forces_hartree_per_bohr = flatten_forces_hartree_per_bohr(force, ucell.nat); } catch (const std::exception& exc) { diff --git a/source/source_relax/test/socket_driver_test.cpp b/source/source_relax/test/socket_driver_test.cpp index 7f8c4255c2a..72478363025 100644 --- a/source/source_relax/test/socket_driver_test.cpp +++ b/source/source_relax/test/socket_driver_test.cpp @@ -22,20 +22,6 @@ #include #include -namespace unitcell -{ -void periodic_boundary_adjustment(Atom* atoms, const ModuleBase::Matrix3& latvec, const int ntype) -{ - for (int it = 0; it < ntype; ++it) - { - for (int ia = 0; ia < atoms[it].na; ++ia) - { - atoms[it].tau[ia] = atoms[it].taud[ia] * latvec; - } - } -} -} // namespace unitcell - namespace { constexpr std::size_t IPI_HEADER_LEN = 12; @@ -244,6 +230,9 @@ class FakeESolver : public ModuleESolver::ESolver void cal_stress(BaseCell&, ModuleBase::matrix& stress) override { stress.create(3, 3); + stress(0, 0) = 2.0; + stress(1, 1) = 2.0; + stress(2, 2) = 2.0; } private: @@ -271,6 +260,7 @@ void initialize_one_atom_cell(UnitCell& ucell) { ucell.lat0 = 1.0; ucell.latvec.Identity(); + ucell.omega = 1.0; ucell.ntype = 1; ucell.nat = 1; ucell.atoms[0].na = 1; From 9b2f3622d291333c237e9592f270d792fc3456cf Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Fri, 14 Aug 2026 03:49:37 +0800 Subject: [PATCH 18/24] fix: publish socket convergence and property metadata --- ...6-08-05-fixed-cell-ipi-hardening-design.md | 4 +- interfaces/ASE_interface/abacuslite/core.py | 27 +++++++++++- source/source_esolver/esolver_ks.cpp | 1 + .../source_relax/test/socket_driver_test.cpp | 42 +++++++++++++++++++ 4 files changed, 72 insertions(+), 2 deletions(-) diff --git a/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md b/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md index 4addb2822d3..cc4a58a8f5e 100644 --- a/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md +++ b/docs/superpowers/specs/2026-08-05-fixed-cell-ipi-hardening-design.md @@ -113,7 +113,9 @@ mistaken transpose cannot pass accidentally. After `runner()`: -- reject a KS frame when SCF did not converge; +- publish a finite KS frame even when SCF did not converge, and set + `scf_converged=false` in the i-PI extras metadata so the external optimizer + or integrator can apply its own acceptance policy; - require finite energy; - require force shape `nat x 3` and finite force entries; - when stress is enabled, require a finite `3 x 3` stress tensor, positive diff --git a/interfaces/ASE_interface/abacuslite/core.py b/interfaces/ASE_interface/abacuslite/core.py index e1cae67fc82..da9c154572b 100644 --- a/interfaces/ASE_interface/abacuslite/core.py +++ b/interfaces/ASE_interface/abacuslite/core.py @@ -697,7 +697,11 @@ def calculate(self, atoms=None, properties=None, system_changes=None): metadata = self._decode_socket_metadata(results.pop('morebytes', None)) self._last_socket_metadata = metadata if metadata is None: - present = set(self._active_properties) + if set(self._active_properties) != {'energy'}: + raise ValueError( + 'ABACUS socket response omitted property metadata; refusing ' + 'to infer forces or stress from fixed-wire padding') + present = {'energy'} self.last_scf_converged = None else: present = set(metadata['present']) @@ -992,6 +996,27 @@ def test_socketio_metadata_marks_padding_absent(self): self.assertEqual(decoded['present'], ('energy',)) self.assertFalse(decoded['scf_converged']) + def test_socketio_legacy_response_cannot_infer_force_from_padding(self): + class LegacyServer: + def calculate(self, atoms): + return { + 'energy': 1.0, + 'forces': np.zeros((len(atoms), 3)), + 'virial': np.zeros((3, 3)), + 'morebytes': b'', + } + + calc = object.__new__(AbacusSocketIO) + calc.variable_cell = False + calc._property_constraints = {} + calc._active_properties = ('energy', 'forces') + calc._reference_cell = None + calc.atoms = None + calc.server = LegacyServer() + atoms = Atoms('Si') + with self.assertRaisesRegex(ValueError, 'refusing to infer forces'): + calc.calculate(atoms=atoms, properties=('forces',), system_changes=()) + def test_socketio_requested_disabled_property_is_rejected(self): calc = object.__new__(AbacusSocketIO) calc._property_constraints = {'forces': False, 'stress': False} diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index f1e7c240be0..a5b1c2fdc49 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -166,6 +166,7 @@ void ESolver_KS::runner(BaseCell& basecell, const int istep) // 7) after scf this->after_scf(ucell, istep, conv_esolver); + this->conv_esolver = conv_esolver; ModuleBase::timer::end(this->classname, "runner"); return; diff --git a/source/source_relax/test/socket_driver_test.cpp b/source/source_relax/test/socket_driver_test.cpp index 72478363025..f3aefe60e3c 100644 --- a/source/source_relax/test/socket_driver_test.cpp +++ b/source/source_relax/test/socket_driver_test.cpp @@ -8,6 +8,7 @@ #include "for_test.h" #include +#include #include #include #include @@ -516,3 +517,44 @@ TEST(SocketDriverTest, EnergyAndStressFrameDoesNotAdvertiseForce) EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"forces\""))); EXPECT_NE(0.0, response.virial_wire_hartree[0]); } + +TEST(SocketDriverTest, EnergyAndForceFrameAdvertisesOnlyForce) +{ + ForceResponse response; + const DriverResult result = run_driver_frame( + true, true, false, + [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); + + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"forces\"]")); + EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"stress\""))); + EXPECT_THAT(response.forces_hartree_per_bohr, + testing::ElementsAre(0.0, 0.0, 0.0)); + EXPECT_THAT(response.virial_wire_hartree, + testing::ElementsAre(0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0)); +} + +TEST(SocketDriverTest, EnergyForceAndStressFrameAdvertisesBothDerivatives) +{ + ForceResponse response; + const DriverResult result = run_driver_frame( + true, true, true, + [&](const int fd) { + send_header(fd, "GETFORCE"); + response = read_force_response(fd); + }); + + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(response.extra, + testing::HasSubstr("\"present\":[\"energy\",\"forces\",\"stress\"]")); + EXPECT_EQ(3u, response.forces_hartree_per_bohr.size()); + EXPECT_TRUE(std::isfinite(response.forces_hartree_per_bohr[0])); + EXPECT_NE(0.0, response.virial_wire_hartree[0]); +} From c3f4b699001236432b0c6eb6cff5d6ff31e7bcce Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Fri, 14 Aug 2026 04:19:31 +0800 Subject: [PATCH 19/24] fix: align fixed socket convergence and property contract --- docs/advanced/input_files/input-main.md | 2 +- docs/advanced/interface/ase.md | 2 ++ docs/parameters.yaml | 2 +- source/source_esolver/esolver_ks_lcao.cpp | 2 +- source/source_relax/test/socket_driver_test.cpp | 6 ++++-- 5 files changed, 9 insertions(+), 5 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 36d7d9d8455..c7e777b70f0 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -604,7 +604,7 @@ - host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. - path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. -- Socket mode always computes energy. Force and stress extraction follow `cal_force` and `cal_stress` independently. A disabled property is sent only as fixed-protocol padding and is marked absent in i-PI extras metadata, so a client must not interpret that padding as a physical zero. If SCF does not converge, ABACUS returns the available frame with `scf_converged=false` metadata; the external optimizer or integrator decides whether to accept it. +- Socket mode always computes energy. Force and stress extraction follow `cal_force` and `cal_stress` independently. A disabled property is sent only as fixed-protocol padding and is marked absent in the ABACUS i-PI extras metadata, so a client must not interpret that padding as a physical zero. The metadata extension is required for safe force/stress handling; a legacy empty-extras response is accepted only for energy-only use, and a generic client that ignores extras cannot distinguish padding from a computed zero. If SCF does not converge, ABACUS returns the available frame with `scf_converged=false` metadata; the external optimizer or integrator decides whether to accept it. - **Default**: False ### esolver_type diff --git a/docs/advanced/interface/ase.md b/docs/advanced/interface/ase.md index 57a3f5841e5..0d753ce6802 100644 --- a/docs/advanced/interface/ase.md +++ b/docs/advanced/interface/ase.md @@ -166,6 +166,8 @@ There are two endpoint styles: Calling `atoms.get_potential_energy()` does not force a force or stress calculation. If a requested property was disabled, ASE raises `PropertyNotImplementedError`; zero-filled i-PI padding is never treated as a physical result. When SCF does not converge, `AbacusSocketIO.last_scf_converged` is set to `False` and the caller decides whether to continue or stop. +The ABACUS metadata extension is required to expose force/stress presence safely. If a legacy client returns an empty extras field, the adapter accepts only an energy-only response and refuses to infer forces or stress from the fixed-wire padding. Generic i-PI/ASE clients that ignore ABACUS extras cannot distinguish mandatory padding from a computed zero; use `AbacusSocketIO` or another metadata-aware client when requesting optional properties. + A socket calculator owns one ABACUS process initialized from one fixed `INPUT`/`STRU` setup. Reuse the same `AbacusSocketIO` instance only for position updates under the same electronic-structure settings and the same cell. Do not change `kpts`, `kspacing`, `nspin`, `basis_type`, `basissets`, pseudopotentials, species, atom count, cell, or other core `INPUT`/`STRU` parameters through an existing socket calculator; create a new `AbacusSocketIO` instance and a new ABACUS client process for those changes. `AbacusSocketIO` rejects cell changes before sending them to ABACUS, and the ABACUS socket driver also checks incoming POSDATA cells against the initial `STRU` cell and exits if they differ. In socket mode, ABACUS keeps one client process alive. All SCF evaluations produced by the same `AbacusSocketIO` instance are appended to the same `OUT.ABACUS/running_scf.log`, because the ABACUS calculation type remains `scf`. The authoritative per-step energy and force results are returned through the i-PI socket to ASE. Use ASE trajectory and optimizer log files, such as `BFGS(atoms, trajectory="opt.traj", logfile="opt.log")`, when each optimizer or MD step should be saved separately. Treat `running_scf.log` mainly as the ABACUS diagnostic log for the socket client, not as one independent FileIO result per structure. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index 138832d008d..113fe3d3eb5 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -49,7 +49,7 @@ parameters: * path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. - Socket mode always computes energy. Force and stress extraction follows cal_force and cal_stress independently; disabled properties are sent as protocol padding and marked absent in the i-PI extras metadata, not reported as physical zero values. A non-converged SCF step is returned with scf_converged=false metadata so an external driver can choose its policy. + Socket mode always computes energy. Force and stress extraction follows cal_force and cal_stress independently; disabled properties are sent as protocol padding and marked absent in the ABACUS i-PI extras metadata, not reported as physical zero values. This metadata extension is required for safe optional-property handling: a legacy response with empty extras is accepted only for energy-only use, while a generic client that ignores extras cannot distinguish padding from a computed zero. A non-converged SCF step is returned with scf_converged=false metadata so an external driver can choose its policy. default_value: "False" unit: "" availability: "" diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 422e2f9fefa..8558c2e14c1 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -566,7 +566,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep, const this->orb_, this->pw_wfc, this->pw_rho, this->pw_big, this->sf, this->pw_rhod, this->locpp.vloc, this->solvent, this->rdmft_solver, this->deepks, this->exx_nao, - this->conv_esolver, this->scf_nmax_flag, istep); + conv_esolver, this->scf_nmax_flag, istep); //! 3) Clean up RA, which is used to serach for adjacent atoms if (!PARAM.inp.cal_force && !PARAM.inp.cal_stress) diff --git a/source/source_relax/test/socket_driver_test.cpp b/source/source_relax/test/socket_driver_test.cpp index f3aefe60e3c..452c632df1e 100644 --- a/source/source_relax/test/socket_driver_test.cpp +++ b/source/source_relax/test/socket_driver_test.cpp @@ -226,6 +226,7 @@ class FakeESolver : public ModuleESolver::ESolver void cal_force(BaseCell& cell, ModuleBase::matrix& force) override { force.create(cell.nat(), 3); + force(0, 0) = 4.0; } void cal_stress(BaseCell&, ModuleBase::matrix& stress) override @@ -533,7 +534,7 @@ TEST(SocketDriverTest, EnergyAndForceFrameAdvertisesOnlyForce) EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"forces\"]")); EXPECT_THAT(response.extra, testing::Not(testing::HasSubstr("\"stress\""))); EXPECT_THAT(response.forces_hartree_per_bohr, - testing::ElementsAre(0.0, 0.0, 0.0)); + testing::ElementsAre(2.0, 0.0, 0.0)); EXPECT_THAT(response.virial_wire_hartree, testing::ElementsAre(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, @@ -555,6 +556,7 @@ TEST(SocketDriverTest, EnergyForceAndStressFrameAdvertisesBothDerivatives) EXPECT_THAT(response.extra, testing::HasSubstr("\"present\":[\"energy\",\"forces\",\"stress\"]")); EXPECT_EQ(3u, response.forces_hartree_per_bohr.size()); - EXPECT_TRUE(std::isfinite(response.forces_hartree_per_bohr[0])); + EXPECT_THAT(response.forces_hartree_per_bohr, + testing::ElementsAre(2.0, 0.0, 0.0)); EXPECT_NE(0.0, response.virial_wire_hartree[0]); } From 76ab2e90706332d0501a381266d5dafe684949ef Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 6 Sep 2026 21:42:19 +0800 Subject: [PATCH 20/24] Fix Makefile linkage for fixed-cell i-PI socket driver --- source/Makefile.Objects | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c6d2ec116e0..29ef4c8c86e 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -535,6 +535,8 @@ OBJS_PW=fft_bundle.o\ pw_op.o\ OBJS_RELAXATION=relax_data.o\ + socket_ipi.o\ + socket_driver.o\ cg_base.o\ bfgs_basic.o\ relax_driver.o\ From 8ff24b26805e61dc14e28ac6847887b2958a313b Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Sun, 6 Sep 2026 23:18:36 +0800 Subject: [PATCH 21/24] Fix CUSOLVERMp CUDA warning include --- source/source_hsolver/kernels/cuda/diag_cusolvermp.cu | 1 + 1 file changed, 1 insertion(+) diff --git a/source/source_hsolver/kernels/cuda/diag_cusolvermp.cu b/source/source_hsolver/kernels/cuda/diag_cusolvermp.cu index 6d184fd59ef..34c97b20508 100644 --- a/source/source_hsolver/kernels/cuda/diag_cusolvermp.cu +++ b/source/source_hsolver/kernels/cuda/diag_cusolvermp.cu @@ -1,6 +1,7 @@ #ifdef __CUSOLVERMP #include "diag_cusolvermp.cuh" #include "source_base/module_device/device_check.h" +#include "source_base/global_function.h" #include From e8e0f9ae6ad741f0e6d3dafe28ad498c3d9b45fa Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Mon, 7 Sep 2026 11:29:30 +0800 Subject: [PATCH 22/24] Sync socket parameter metadata with generated docs --- docs/advanced/input_files/input-main.md | 7 ++++--- source/source_io/module_parameter/read_inp_sys.cpp | 9 +++++++-- 2 files changed, 11 insertions(+), 5 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e14a7e43bf5..90bf08dfca8 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -648,7 +648,8 @@ - host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. - path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument. -- Socket mode always computes energy. Force and stress extraction follow `cal_force` and `cal_stress` independently. A disabled property is sent only as fixed-protocol padding and is marked absent in the ABACUS i-PI extras metadata, so a client must not interpret that padding as a physical zero. The metadata extension is required for safe force/stress handling; a legacy empty-extras response is accepted only for energy-only use, and a generic client that ignores extras cannot distinguish padding from a computed zero. If SCF does not converge, ABACUS returns the available frame with `scf_converged=false` metadata; the external optimizer or integrator decides whether to accept it. + + Socket mode always computes energy. Force and stress extraction follows cal_force and cal_stress independently; disabled properties are sent as protocol padding and marked absent in the ABACUS i-PI extras metadata, not reported as physical zero values. This metadata extension is required for safe optional-property handling: a legacy response with empty extras is accepted only for energy-only use, while a generic client that ignores extras cannot distinguish padding from a computed zero. A non-converged SCF step is returned with scf_converged=false metadata so an external driver can choose its policy. - **Default**: False ### esolver_type @@ -700,7 +701,7 @@ - **Type**: Boolean - **Description**: If set to True, calculate the force at the end of the electronic iteration. -- In `socket_driver` mode, this controls whether forces are calculated and advertised; it is not forced on by the socket protocol. + In socket_driver mode, this flag controls whether the returned frame advertises forces; it is not forced on by the socket protocol. - **Default**: False ### kpar @@ -816,7 +817,7 @@ - **Type**: Boolean - **Description**: If set to True, calculate the stress at the end of the electronic iteration. -- In `socket_driver` mode, this independently controls whether stress/virial is calculated and advertised. + In socket_driver mode, this flag independently controls whether the returned frame advertises stress/virial. - **Default**: False ### diago_proc diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index b93df595b6e..795a6208c99 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -179,6 +179,9 @@ void ReadInput::item_system() * host:port, for example localhost:31415 or 127.0.0.1:31415, opens a TCP connection to that host and port. Use this when the i-PI server listens on a TCP port. * path:UNIX, for example /tmp/ipi_abacus_si:UNIX, opens a Unix-domain socket at the given filesystem path. The :UNIX suffix tells ABACUS that the preceding value is a local socket path rather than a TCP host name. This form only works on the same machine. When using the ASE AbacusSocketIO interface, this environment variable is set automatically from the port or unixsocket calculator argument.)"; + item.description += R"( + +Socket mode always computes energy. Force and stress extraction follows cal_force and cal_stress independently; disabled properties are sent as protocol padding and marked absent in the ABACUS i-PI extras metadata, not reported as physical zero values. This metadata extension is required for safe optional-property handling: a legacy response with empty extras is accepted only for energy-only use, while a generic client that ignores extras cannot distinguish padding from a computed zero. A non-converged SCF step is returned with scf_converged=false metadata so an external driver can choose its policy.)"; item.default_value = "False"; read_sync_bool(input.socket_driver); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -323,7 +326,8 @@ When using the ASE AbacusSocketIO interface, this environment variable is set au item.annotation = "if calculate the force at the end of the electronic iteration"; item.category = "System variables"; item.type = "Boolean"; - item.description = "If set to True, calculate the force at the end of the electronic iteration."; + item.description = R"(If set to True, calculate the force at the end of the electronic iteration. +In socket_driver mode, this flag controls whether the returned frame advertises forces; it is not forced on by the socket protocol.)"; item.default_value = "False"; item.reset_value = [](const Input_Item& item, Parameter& para) { std::vector use_force = {"cell-relax", "relax", "md"}; @@ -627,7 +631,8 @@ For `basis_type=lcao_in_pw`, `init_wfc` is automatically set to `nao`. item.annotation = "calculate the stress or not"; item.category = "System variables"; item.type = "Boolean"; - item.description = "If set to True, calculate the stress at the end of the electronic iteration."; + item.description = R"(If set to True, calculate the stress at the end of the electronic iteration. +In socket_driver mode, this flag independently controls whether the returned frame advertises stress/virial.)"; item.default_value = "False"; item.reset_value = [](const Input_Item& item, Parameter& para) { if (para.input.calculation == "md") From 29b18ac1d743b73f9061433fa7878c4b4f538c5a Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Mon, 7 Sep 2026 14:37:48 +0800 Subject: [PATCH 23/24] Docs: clarify abacuslite README scope --- interfaces/ASE_interface/README.md | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/interfaces/ASE_interface/README.md b/interfaces/ASE_interface/README.md index a76e68e177e..3b1567364bf 100644 --- a/interfaces/ASE_interface/README.md +++ b/interfaces/ASE_interface/README.md @@ -7,15 +7,17 @@ abacuslite is a lightweight plugin for ABACUS (Atomic-orbital Based Ab-initio Co ### Key Features - **Lightweight Design**: Implemented as a plugin, no need to modify ASE core code -- **Version Compatibility**: No longer restricted to specific ASE versions, works with most ASE versions +- **Version Compatibility**: Supports ASE versions satisfying the package requirement `ase>=3.22` - **ASE Integration**: Uses ASE as the running platform, making ABACUS a callable calculator within it -- **Function Support**: Currently only supports SCF (Self-Consistent Field) functionality, returning energy, forces, stress, etc. +- **Function Support**: Provides SCF-based energy, force, and stress evaluations through ASE. ASE can use these evaluations for relaxation, molecular dynamics, NEB, band-structure, and density-of-states workflows. +- **Socket Support**: `AbacusSocketIO` provides fixed-cell i-PI socket calculations, with energy always available and forces/stress enabled independently when requested. ## Installation -Installation is very simple, just execute the following command in the project root directory: +Install the plugin from the ASE interface directory: ```bash +cd interfaces/ASE_interface pip install . ``` @@ -32,9 +34,9 @@ Please refer to the example scripts in the `examples` folder. Recommended learni 7. **constraintmd.py** - Constrained molecular dynamics simulation 8. **metadynamics.py** - Metadynamics simulation 9. **neb.py** - Nudged Elastic Band (NEB) calculation -10. **socketio.py** - ASE optimization with the AbacusSocketIO calculator running ABACUS as an i-PI socket client +10. **socketio.py** - Fixed-cell ASE optimization with `AbacusSocketIO`, running ABACUS as an i-PI socket client -More usage examples will be provided in future versions. +The regular `Abacus` calculator runs one ABACUS calculation for each ASE property evaluation. ASE controls the relaxation, molecular-dynamics, and other workflow steps. The socket calculator reuses one ABACUS process for position updates, while the cell and electronic-structure settings remain fixed for that calculator instance. ## Authors @@ -49,10 +51,10 @@ Thanks to the ABACUS development team for their support and contributions. ## License -[Fill in according to the actual project license] +This interface is distributed under the [GNU Lesser General Public License v3.0 or later](../../LICENSE). ## Contact If you have any questions or suggestions, please contact us through: -- GitHub: [deepmodeling/abacus-develop](https://github.com/deepmodeling/abacus-develop) \ No newline at end of file +- GitHub: [deepmodeling/abacus-develop](https://github.com/deepmodeling/abacus-develop) From 46ec7447c2a49264345e078989576c2c4ce3cd96 Mon Sep 17 00:00:00 2001 From: jianrui geng Date: Mon, 7 Sep 2026 14:41:19 +0800 Subject: [PATCH 24/24] Docs: clarify abacuslite license and examples --- interfaces/ASE_interface/README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/interfaces/ASE_interface/README.md b/interfaces/ASE_interface/README.md index 3b1567364bf..543d4774836 100644 --- a/interfaces/ASE_interface/README.md +++ b/interfaces/ASE_interface/README.md @@ -34,7 +34,8 @@ Please refer to the example scripts in the `examples` folder. Recommended learni 7. **constraintmd.py** - Constrained molecular dynamics simulation 8. **metadynamics.py** - Metadynamics simulation 9. **neb.py** - Nudged Elastic Band (NEB) calculation -10. **socketio.py** - Fixed-cell ASE optimization with `AbacusSocketIO`, running ABACUS as an i-PI socket client +10. **soc.py** - Noncollinear spin-orbit coupling calculation +11. **socketio.py** - Fixed-cell ASE optimization with `AbacusSocketIO`, running ABACUS as an i-PI socket client The regular `Abacus` calculator runs one ABACUS calculation for each ASE property evaluation. ASE controls the relaxation, molecular-dynamics, and other workflow steps. The socket calculator reuses one ABACUS process for position updates, while the cell and electronic-structure settings remain fixed for that calculator instance. @@ -51,7 +52,7 @@ Thanks to the ABACUS development team for their support and contributions. ## License -This interface is distributed under the [GNU Lesser General Public License v3.0 or later](../../LICENSE). +The applicable license terms are provided in the repository [LICENSE](../../LICENSE). ## Contact