diff --git a/source/Makefile.Objects b/source/Makefile.Objects index fc1306b5c07..9e68b3674b6 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -537,7 +537,10 @@ OBJS_PW=fft_bundle.o\ OBJS_RELAXATION=relax_data.o\ socket_ipi.o\ socket_frame.o\ + socket_frame_utils.o\ socket_driver.o\ + socket_utils.o\ + socket_handlers.o\ cg_base.o\ bfgs_basic.o\ relax_driver.o\ diff --git a/source/source_base/parallel_reduce.cpp b/source/source_base/parallel_reduce.cpp index eaafd7bf6ca..62ac884f0a2 100644 --- a/source/source_base/parallel_reduce.cpp +++ b/source/source_base/parallel_reduce.cpp @@ -112,6 +112,7 @@ template void Parallel_Reduce::reduce_min(int&); template void Parallel_Reduce::reduce_min(float&); template void Parallel_Reduce::reduce_min(double&); +template void Parallel_Reduce::reduce_max(int&); template void Parallel_Reduce::reduce_max(float&); template void Parallel_Reduce::reduce_max(double&); diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index b32dda901cc..691f91b20d0 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -3,7 +3,10 @@ add_library( OBJECT relax_data.cpp socket_ipi.cpp + socket_frame_utils.cpp socket_frame.cpp + socket_utils.cpp + socket_handlers.cpp socket_driver.cpp cg_base.cpp relax_driver.cpp diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index a5bf94cde21..d63c83441d5 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -1,394 +1,53 @@ #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" +#include "source_relax/socket_handlers.h" +#include "source_relax/socket_utils.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 -#include -#include -#include #include -#include -#include -#include +#include #include -#include -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() -{ -#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) -{ -#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_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 - 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()); - bcast_socket_int(size); - if (!is_root()) - { - value.resize(static_cast(size)); - } - if (size > 0) - { - bcast_socket_chars(&value[0], size); - } -} - -void quit_if_root_io_failed(int root_failed, std::string root_message) -{ - bcast_socket_int(root_failed); - bcast_socket_string(root_message); - if (root_failed != 0) - { - ModuleBase::WARNING_QUIT("ABACUS socket", root_message.empty() ? "i-PI socket I/O failed" : root_message); - } -} - -std::string bcast_header(std::string header) -{ - bcast_socket_string(header); - return header; -} - -std::string socket_address() -{ - const char* env = std::getenv("ABACUS_SOCKET_ADDRESS"); - if (env == nullptr || std::string(env).empty()) - { - return "localhost:31415"; - } - return std::string(env); -} +using SocketUtils::ComputedFrame; +using SocketUtils::DriverState; +using SocketUtils::ipi_cell_bohr_from_unitcell; +using SocketUtils::is_root; +using SocketUtils::quit_if_root_failed; +using SocketUtils::address; +using SocketHandlers::DriverContext; -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_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()) - { - 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; -} - -double unchanged_cell_tolerance(const SocketFrame::Matrix9& cell) +namespace { - 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) +void connect_on_root(IpiSocket& socket, std::ofstream& ofs_running) { - if (positions_bohr.size() != static_cast(3 * ucell.nat)) - { - ModuleBase::WARNING_QUIT("ABACUS socket", "POSDATA atom count does not match STRU."); - } - - int iat = 0; - for (int it = 0; it < ucell.ntype; ++it) + int io_failed = 0; + std::string io_message; + if (is_root()) { - Atom* atom = &ucell.atoms[it]; - for (int ia = 0; ia < atom->na; ++ia) + try { - 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; + const std::string endpoint = address(); + ofs_running << " ABACUS socket driver connecting to i-PI endpoint " << endpoint << std::endl; + socket.connect(endpoint); } - } - 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, const int nat) -{ - if (nat < 0 || force.nr != nat || 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) + catch (const std::exception& exc) { - 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; + io_failed = 1; + io_message = exc.what(); } } - return out; + quit_if_root_failed(io_failed, io_message); } -SocketFrame::Matrix9 matrix9_from_stress(const ModuleBase::matrix& stress) +void log_peer_closed(std::ofstream& ofs_running) { - 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) + if (is_root()) { - for (int column = 0; column < 3; ++column) - { - values[3 * row + column] = stress(row, column); - } + ofs_running << " ABACUS socket driver exiting after peer closed connection" << std::endl; } - return values; -} - -std::vector vector_from_matrix9(const SocketFrame::Matrix9& values) -{ - return std::vector(values.begin(), values.end()); } } // namespace @@ -408,456 +67,50 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, try { - int io_failed = 0; - std::string io_message; - if (is_root()) - { - try - { - const std::string address = socket_address(); - 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(); - } - } - quit_if_root_io_failed(io_failed, io_message); - - DriverState state = DriverState::NeedInit; - int istep = 0; - const int nat_return = ucell.nat; - ComputedFrame published; + connect_on_root(socket, ofs_running); - const std::vector reference_cell = ipi_cell_bohr_from_unitcell(ucell); - bool checked_initial_positions = false; + DriverContext context; + context.esolver = p_esolver; + context.ucell = &ucell; + context.inp = &inp; + context.state = DriverState::NeedInit; + context.nat_return = ucell.nat; + context.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 IpiSocketClosed&) - { - 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) - { - io_failed = 1; - io_message = exc.what(); - } - } - quit_if_root_io_failed(io_failed, io_message); - header = bcast_header(header); - + const std::string header + = SocketHandlers::read_header_bcast(socket, context.state); if (header.empty()) { - if (is_root()) - { - ofs_running << " ABACUS socket driver exiting after peer closed connection" << std::endl; - } + log_peer_closed(ofs_running); break; } else if (header == "STATUS") { - io_failed = 0; - io_message.clear(); - if (is_root()) - { - try - { - if (state == DriverState::HasData) - { - socket.write_header("HAVEDATA"); - } - else if (state == DriverState::Ready) - { - socket.write_header("READY"); - } - else - { - socket.write_header("NEEDINIT"); - } - } - catch (const std::exception& exc) - { - io_failed = 1; - io_message = exc.what(); - } - } - quit_if_root_io_failed(io_failed, io_message); + SocketHandlers::handle_status(socket, context.state); } else if (header == "INIT") { - std::int32_t rid = 0; - std::int32_t nbytes = 0; - std::string params; - io_failed = 0; - io_message.clear(); - if (is_root()) - { - if (state != DriverState::NeedInit) - { - io_failed = 1; - io_message = "INIT requires NEEDINIT state"; - } - else - { - try - { - 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)); - } - } - catch (const std::exception& exc) - { - io_failed = 1; - io_message = exc.what(); - } - } - } - quit_if_root_io_failed(io_failed, io_message); - bcast_socket_int32(rid); - bcast_socket_int32(nbytes); - if (nbytes > 0 && is_root()) - { - ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; - } - state = DriverState::Ready; - if (is_root()) - { - ofs_running << " ABACUS socket INIT replica " << rid << std::endl; - } + SocketHandlers::handle_init(socket, context, ofs_running); } else if (header == "POSDATA") { - 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()) - { - if (state != DriverState::Ready) - { - io_failed = 1; - io_message = "POSDATA requires READY state"; - } - else - { - try - { - 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; - } - } - } - catch (const std::exception& exc) - { - io_failed = 1; - io_message = exc.what(); - } - } - } - quit_if_root_io_failed(io_failed, io_message); - 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()) - { - 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); - - 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."); - } - 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."); - } - } - - 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, " - << computed.energy_hartree << " Ha" << std::endl; - } - ModuleBase::matrix force; - if (inp.cal_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, ucell.nat); - } - 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; - state = DriverState::HasData; + SocketHandlers::handle_posdata(socket, context, ofs_running); } else if (header == "GETFORCE") { - io_failed = 0; - io_message.clear(); - if (is_root()) - { - 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(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) - { - io_failed = 1; - io_message = exc.what(); - } - } - quit_if_root_io_failed(io_failed, io_message); - published = ComputedFrame(); - state = DriverState::Ready; + SocketHandlers::handle_getforce(socket, context); } else if (header == "EXIT") { - if (is_root()) - { - ofs_running << " ABACUS socket driver received i-PI EXIT" << std::endl; - } + SocketHandlers::handle_exit(ofs_running); break; } else { - if (is_root()) - { - io_failed = 1; - io_message = "unknown i-PI header: " + header; - } - quit_if_root_io_failed(io_failed, io_message); + quit_if_root_failed(is_root() ? 1 : 0, + is_root() ? "unknown i-PI header: " + header : ""); } } } diff --git a/source/source_relax/socket_frame.cpp b/source/source_relax/socket_frame.cpp index 9f125310761..1a009879ddb 100644 --- a/source/source_relax/socket_frame.cpp +++ b/source/source_relax/socket_frame.cpp @@ -1,163 +1,202 @@ #include "socket_frame.h" +#include "source_relax/socket_frame_utils.h" + #include #include #include +using namespace FrameUtils; + +namespace SocketFrame +{ namespace { -const int MATRIX_DIMENSION = 3; -const int MAX_JACOBI_SWEEPS = 32; +CellValidation make_failed_cell_validation() +{ + 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); + return result; +} -bool is_finite_matrix(const SocketFrame::Matrix9& values) +bool validate_cell_entries(const Matrix9& cell_wire, + const Matrix9& inverse_wire, + std::string& message) { - for (std::size_t index = 0; index < values.size(); ++index) + if (!is_finite_matrix(cell_wire) || !is_finite_matrix(inverse_wire)) { - if (!std::isfinite(values[index])) - { - return false; - } + message = "cell and received inverse entries must be finite"; + return false; } return true; } -double column_norm_squared(const SocketFrame::Matrix9& values, int column) +bool validate_cell_tolerances(const double max_condition_number, + const double inverse_absolute_tolerance, + const double inverse_relative_tolerance, + std::string& message) { - double norm_squared = 0.0; - for (int row = 0; row < MATRIX_DIMENSION; ++row) + 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) { - const double value = values[row * MATRIX_DIMENSION + column]; - norm_squared += value * value; + message = "cell validation tolerances must be finite and nonnegative"; + return false; } - return norm_squared; + return true; } -double column_dot(const SocketFrame::Matrix9& values, int first, int second) +bool compute_cell_scale(const Matrix9& cell_wire, double& scale, Matrix9& scaled_cell, + std::string& message) { - double dot = 0.0; - for (int row = 0; row < MATRIX_DIMENSION; ++row) + 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) + { + message = "cell determinant must be positive"; + return false; + } + for (std::size_t index = 0; index < cell_wire.size(); ++index) { - dot += values[row * MATRIX_DIMENSION + first] * values[row * MATRIX_DIMENSION + second]; + scaled_cell[index] = cell_wire[index] / scale; } - return dot; + return true; } -bool columns_are_orthogonal(const SocketFrame::Matrix9& values) +bool compute_cell_determinant(const double scale, + const Matrix9& scaled_cell, + double& determinant_bohr3, + std::string& message) { - 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; - } + const long double determinant_scaled = scaled_determinant(scaled_cell); + if (determinant_scaled <= 0.0L) + { + message = "cell determinant must be positive"; + return false; + } + 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())) + { + message = "cell determinant is not representable as a finite double"; + return false; + } + determinant_bohr3 = static_cast(determinant); + if (!std::isfinite(determinant_bohr3) || determinant_bohr3 <= 0.0) + { + message = "cell determinant is not representable as a positive finite double"; + return false; } return true; } -void rotate_columns(SocketFrame::Matrix9& values, int first, int second, double cosine, double sine) +bool compute_cell_svd(const Matrix9& scaled_cell, + double singular_values[kMatrixDimension], + Matrix9& orthogonal_columns, + Matrix9& right_vectors, + std::string& message) { - for (int row = 0; row < MATRIX_DIMENSION; ++row) + orthogonal_columns = scaled_cell; + if (!one_sided_jacobi(orthogonal_columns, right_vectors)) { - 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; + message = "cell singular-value iteration did not converge"; + return false; } + for (int column = 0; column < kMatrixDimension; ++column) + { + singular_values[column] = std::sqrt(column_norm_squared(orthogonal_columns, column)); + } + return true; } -bool one_sided_jacobi(SocketFrame::Matrix9& columns, SocketFrame::Matrix9& right_vectors) +bool compute_condition_number(const double singular_values[kMatrixDimension], + const double max_condition_number, + double& condition_number_2, + std::string& message) { - 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}}; + double largest_singular = 0.0; + double smallest_singular = std::numeric_limits::infinity(); + for (int column = 0; column < kMatrixDimension; ++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)) + { + message = "cell is singular"; + return false; + } + condition_number_2 = largest_singular / smallest_singular; + if (!std::isfinite(condition_number_2) + || condition_number_2 >= max_condition_number) + { + message = "cell condition number is not below the configured maximum"; + return false; + } + return true; +} - for (int sweep = 0; sweep < MAX_JACOBI_SWEEPS; ++sweep) +void compute_cell_inverse(const double scale, + const double singular_values[kMatrixDimension], + const Matrix9& orthogonal_columns, + const Matrix9& right_vectors, + Matrix9& computed_inverse) +{ + for (int row = 0; row < kMatrixDimension; ++row) { - for (int pair = 0; pair < 3; ++pair) + for (int column = 0; column < kMatrixDimension; ++column) { - 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) + long double inverse_value = 0.0L; + for (int singular = 0; singular < kMatrixDimension; ++singular) { - continue; + const long double sigma = singular_values[singular]; + inverse_value + += static_cast(right_vectors[row * kMatrixDimension + singular]) + * orthogonal_columns[column * kMatrixDimension + singular] + / (static_cast(scale) * sigma * sigma); } - - 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; + computed_inverse[row * kMatrixDimension + column] + = static_cast(inverse_value); } } - 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, - bool transpose_inverse) +bool check_inverse_consistency(const Matrix9& cell_wire, + const Matrix9& inverse_wire, + const double condition_number_2, + const double inverse_absolute_tolerance, + const double inverse_relative_tolerance, + double& inverse_residual, + std::string& message) { - long double maximum = 0.0L; - for (int row = 0; row < MATRIX_DIMENSION; ++row) + const double direct_inverse_residual + = received_inverse_residual(cell_wire, inverse_wire, false); + const double transposed_inverse_residual + = received_inverse_residual(cell_wire, inverse_wire, true); + inverse_residual = std::min(direct_inverse_residual, transposed_inverse_residual); + const double residual_limit + = inverse_absolute_tolerance + + inverse_relative_tolerance * condition_number_2 + * std::numeric_limits::epsilon(); + if (!std::isfinite(inverse_residual) || inverse_residual > residual_limit) { - for (int column = 0; column < MATRIX_DIMENSION; ++column) - { - long double product = 0.0L; - for (int inner = 0; inner < MATRIX_DIMENSION; ++inner) - { - const int inverse_index = transpose_inverse - ? column * MATRIX_DIMENSION + inner - : inner * MATRIX_DIMENSION + column; - product += static_cast(cell[row * MATRIX_DIMENSION + inner]) - * inverse[inverse_index]; - } - const long double expected = row == column ? 1.0L : 0.0L; - maximum = std::max(maximum, std::fabs(product - expected)); - } + message = "received cell inverse is inconsistent with the cell"; + return false; } - return static_cast(maximum); + return true; } } // namespace -namespace SocketFrame -{ Matrix9 transpose_matrix9(const Matrix9& values) { return {{values[0], values[3], values[6], @@ -171,125 +210,51 @@ CellValidation validate_ipi_cell(const Matrix9& cell_wire, 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); + CellValidation result = make_failed_cell_validation(); - if (!is_finite_matrix(cell_wire) || !is_finite_matrix(inverse_wire)) + if (!validate_cell_entries(cell_wire, inverse_wire, result.message)) { - 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) + if (!validate_cell_tolerances(max_condition_number, + inverse_absolute_tolerance, + inverse_relative_tolerance, + result.message)) { - 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) + if (!compute_cell_scale(cell_wire, scale, scaled_cell, result.message)) { - 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())) + if (!compute_cell_determinant(scale, scaled_cell, result.determinant_bohr3, result.message)) { - 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; + double singular_values[kMatrixDimension]; + Matrix9 orthogonal_columns; 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)) + if (!compute_cell_svd(scaled_cell, singular_values, orthogonal_columns, right_vectors, + result.message)) { - 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) + if (!compute_condition_number(singular_values, max_condition_number, + result.condition_number_2, result.message)) { - result.message = "cell condition number is not below 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); - } - } - - const double direct_inverse_residual - = received_inverse_residual(cell_wire, inverse_wire, false); - const double transposed_inverse_residual - = received_inverse_residual(cell_wire, inverse_wire, true); - result.inverse_residual = std::min(direct_inverse_residual, transposed_inverse_residual); - 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) + compute_cell_inverse(scale, singular_values, orthogonal_columns, right_vectors, + result.computed_inverse_wire_bohr_inv); + if (!check_inverse_consistency(cell_wire, inverse_wire, result.condition_number_2, + inverse_absolute_tolerance, inverse_relative_tolerance, + result.inverse_residual, result.message)) { - result.message = "received cell inverse is inconsistent with the cell"; return result; } @@ -344,82 +309,139 @@ bool checked_position_count(std::int32_t nat_socket, return true; } -VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, - double volume_bohr3, - double antisymmetric_absolute_tolerance, - double antisymmetric_relative_tolerance) +namespace +{ +VirialConversion make_failed_virial_conversion() { VirialConversion result; result.ok = false; result.message.clear(); result.wire_virial_hartree.fill(0.0); result.max_antisymmetric_component = 0.0; + return result; +} +bool validate_virial_inputs(const Matrix9& stress_ry_per_bohr3, + const double volume_bohr3, + const double antisymmetric_absolute_tolerance, + const double antisymmetric_relative_tolerance, + std::string& message) +{ if (!is_finite_matrix(stress_ry_per_bohr3)) { - result.message = "stress entries must be finite"; - return result; + message = "stress entries must be finite"; + return false; } if (!std::isfinite(volume_bohr3) || volume_bohr3 <= 0.0) { - result.message = "cell volume must be finite and positive"; - return result; + message = "cell volume must be finite and positive"; + return false; } 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; + message = "stress symmetry tolerances must be finite and nonnegative"; + return false; } + return true; +} +bool check_stress_symmetry(const Matrix9& stress_ry_per_bohr3, + const double antisymmetric_absolute_tolerance, + const double antisymmetric_relative_tolerance, + double& max_antisymmetric_component, + std::string& message) +{ 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 row = 0; row < kMatrixDimension; ++row) { - for (int column = row + 1; column < MATRIX_DIMENSION; ++column) + for (int column = row + 1; column < kMatrixDimension; ++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); + = std::fabs(stress_ry_per_bohr3[row * kMatrixDimension + column] + - stress_ry_per_bohr3[column * kMatrixDimension + row]); + max_antisymmetric_component + = std::max(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) + if (!std::isfinite(max_antisymmetric_component) + || max_antisymmetric_component > symmetry_limit) { - result.message = "stress tensor is not symmetric within tolerance"; - return result; + message = "stress tensor is not symmetric within tolerance"; + return false; } + return true; +} +bool compute_symmetric_virial(const Matrix9& stress_ry_per_bohr3, + const double volume_bohr3, + Matrix9& wire_virial_hartree, + std::string& message) +{ Matrix9 virial; - for (int row = 0; row < MATRIX_DIMENSION; ++row) + for (int row = 0; row < kMatrixDimension; ++row) { - for (int column = 0; column < MATRIX_DIMENSION; ++column) + for (int column = 0; column < kMatrixDimension; ++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]); + * (static_cast(stress_ry_per_bohr3[row * kMatrixDimension + column]) + + stress_ry_per_bohr3[column * kMatrixDimension + 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; + message = "converted virial is not representable as finite doubles"; + return false; } - virial[row * MATRIX_DIMENSION + column] = static_cast(converted); + virial[row * kMatrixDimension + column] = static_cast(converted); } } - result.wire_virial_hartree = transpose_matrix9(virial); + wire_virial_hartree = transpose_matrix9(virial); + return true; +} +} // namespace + +VirialConversion make_ipi_virial(const Matrix9& stress_ry_per_bohr3, + double volume_bohr3, + double antisymmetric_absolute_tolerance, + double antisymmetric_relative_tolerance) +{ + VirialConversion result = make_failed_virial_conversion(); + + if (!validate_virial_inputs(stress_ry_per_bohr3, + volume_bohr3, + antisymmetric_absolute_tolerance, + antisymmetric_relative_tolerance, + result.message)) + { + return result; + } + if (!check_stress_symmetry(stress_ry_per_bohr3, + antisymmetric_absolute_tolerance, + antisymmetric_relative_tolerance, + result.max_antisymmetric_component, + result.message)) + { + return result; + } + if (!compute_symmetric_virial(stress_ry_per_bohr3, + volume_bohr3, + result.wire_virial_hartree, + result.message)) + { + return result; + } + result.ok = true; return result; } diff --git a/source/source_relax/socket_frame_utils.cpp b/source/source_relax/socket_frame_utils.cpp new file mode 100644 index 00000000000..d4507572c2e --- /dev/null +++ b/source/source_relax/socket_frame_utils.cpp @@ -0,0 +1,156 @@ +#include "socket_frame_utils.h" + +#include +#include +#include + +namespace FrameUtils +{ +const int kMaxJacobiSweeps = 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 < kMatrixDimension; ++row) + { + const double value = values[row * kMatrixDimension + 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 < kMatrixDimension; ++row) + { + dot += values[row * kMatrixDimension + first] * values[row * kMatrixDimension + 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 < kMatrixDimension; ++row) + { + const int first_index = row * kMatrixDimension + first; + const int second_index = row * kMatrixDimension + 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 < kMaxJacobiSweeps; ++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, + bool transpose_inverse) +{ + long double maximum = 0.0L; + for (int row = 0; row < kMatrixDimension; ++row) + { + for (int column = 0; column < kMatrixDimension; ++column) + { + long double product = 0.0L; + for (int inner = 0; inner < kMatrixDimension; ++inner) + { + const int inverse_index = transpose_inverse + ? column * kMatrixDimension + inner + : inner * kMatrixDimension + column; + product += static_cast(cell[row * kMatrixDimension + inner]) + * inverse[inverse_index]; + } + const long double expected = row == column ? 1.0L : 0.0L; + maximum = std::max(maximum, std::fabs(product - expected)); + } + } + return static_cast(maximum); +} +} // namespace FrameUtils diff --git a/source/source_relax/socket_frame_utils.h b/source/source_relax/socket_frame_utils.h new file mode 100644 index 00000000000..21e22d872f3 --- /dev/null +++ b/source/source_relax/socket_frame_utils.h @@ -0,0 +1,23 @@ +#ifndef ABACUS_SOURCE_RELAX_SOCKET_FRAME_UTILS_H +#define ABACUS_SOURCE_RELAX_SOCKET_FRAME_UTILS_H + +#include "source_relax/socket_frame.h" + +namespace FrameUtils +{ +constexpr int kMatrixDimension = 3; +extern const int kMaxJacobiSweeps; + +bool is_finite_matrix(const SocketFrame::Matrix9& values); +double column_norm_squared(const SocketFrame::Matrix9& values, int column); +double column_dot(const SocketFrame::Matrix9& values, int first, int second); +bool columns_are_orthogonal(const SocketFrame::Matrix9& values); +void rotate_columns(SocketFrame::Matrix9& values, int first, int second, double cosine, double sine); +bool one_sided_jacobi(SocketFrame::Matrix9& columns, SocketFrame::Matrix9& right_vectors); +long double scaled_determinant(const SocketFrame::Matrix9& values); +double received_inverse_residual(const SocketFrame::Matrix9& cell, + const SocketFrame::Matrix9& inverse, + bool transpose_inverse); +} // namespace FrameUtils + +#endif diff --git a/source/source_relax/socket_handlers.cpp b/source/source_relax/socket_handlers.cpp new file mode 100644 index 00000000000..dc333ced087 --- /dev/null +++ b/source/source_relax/socket_handlers.cpp @@ -0,0 +1,518 @@ +#include "socket_handlers.h" + +#include "source_relax/socket_frame.h" +#include "source_base/global_function.h" +#include "source_cell/unitcell.h" +#include "source_esolver/esolver.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace SocketHandlers +{ +using SocketUtils::ComputedFrame; +using SocketUtils::DriverState; +using SocketUtils::kInverseAbsoluteTolerance; +using SocketUtils::kInverseRelativeTolerance; +using SocketUtils::kMaxCellCondition; +using SocketUtils::kMaxInitBytes; +using SocketUtils::kRyToHartree; +using SocketUtils::kStressAbsoluteTolerance; +using SocketUtils::kStressRelativeTolerance; +using SocketUtils::all_ranks_converged; +using SocketUtils::bcast_double_vector; +using SocketUtils::bcast_header; +using SocketUtils::bcast_int32; +using SocketUtils::fail_during_collective_stage; +using SocketUtils::flatten_forces_hartree_per_bohr; +using SocketUtils::ipi_cell_bohr_from_unitcell; +using SocketUtils::is_root; +using SocketUtils::matrix9_from_stress; +using SocketUtils::max_abs_delta; +using SocketUtils::max_wrapped_direct_delta_from_unitcell; +using SocketUtils::properties_extra; +using SocketUtils::quit_if_root_failed; +using SocketUtils::set_positions_from_ipi_bohr; +using SocketUtils::throw_if_any_rank_failed; +using SocketUtils::unchanged_cell_tolerance; +using SocketUtils::vector_from_matrix9; + +std::string read_header_bcast(IpiSocket& socket, const DriverState state) +{ + std::string header; + int io_failed = 0; + std::string io_message; + if (is_root()) + { + try + { + header = socket.read_header(); + } + catch (const IpiSocketClosed&) + { + if (state == DriverState::HasData) + { + io_failed = 1; + io_message = "i-PI peer closed while a computed frame was pending"; + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + quit_if_root_failed(io_failed, io_message); + return bcast_header(header); +} + +void handle_status(IpiSocket& socket, const DriverState state) +{ + int io_failed = 0; + std::string io_message; + if (is_root()) + { + try + { + if (state == DriverState::HasData) + { + socket.write_header("HAVEDATA"); + } + else if (state == DriverState::Ready) + { + socket.write_header("READY"); + } + else + { + socket.write_header("NEEDINIT"); + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + quit_if_root_failed(io_failed, io_message); +} + +void handle_init(IpiSocket& socket, + DriverContext& context, + std::ofstream& ofs_running) +{ + std::int32_t rid = 0; + std::int32_t nbytes = 0; + std::string params; + int io_failed = 0; + std::string io_message; + if (is_root()) + { + if (context.state != DriverState::NeedInit) + { + io_failed = 1; + io_message = "INIT requires NEEDINIT state"; + } + else + { + try + { + 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 > kMaxInitBytes) + { + io_failed = 1; + io_message = "INIT payload exceeds the 1 MiB socket limit"; + } + else if (nbytes > 0) + { + params = socket.read_string(static_cast(nbytes)); + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + } + quit_if_root_failed(io_failed, io_message); + bcast_int32(rid); + bcast_int32(nbytes); + if (nbytes > 0 && is_root()) + { + ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; + } + context.state = DriverState::Ready; + if (is_root()) + { + ofs_running << " ABACUS socket INIT replica " << rid << std::endl; + } +} + +namespace +{ +struct PosdataPayload +{ + SocketFrame::Matrix9 cell = {{0.0}}; + SocketFrame::Matrix9 inv_cell = {{0.0}}; + std::int32_t nat_socket = 0; + std::vector positions; +}; + +// Root rank reads and validates the POSDATA frame; the results are then +// broadcast to all ranks. Calls WARNING_QUIT on protocol/validation failure. +// The READY-state check must stay inside this function so that every rank +// passes through the single quit_if_root_failed collective sequence below; +// a root-only early quit outside would deadlock if any collective were +// added ahead of it. +PosdataPayload read_posdata(IpiSocket& socket, const UnitCell& ucell, const DriverState state) +{ + PosdataPayload payload; + int io_failed = 0; + std::string io_message; + if (is_root()) + { + if (state != DriverState::Ready) + { + io_failed = 1; + io_message = "POSDATA requires READY state"; + } + else + { + try + { + 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(), payload.cell.begin()); + std::copy(inverse_values.begin(), inverse_values.end(), payload.inv_cell.begin()); + payload.nat_socket = socket.read_int32(); + const SocketFrame::CellValidation validation + = SocketFrame::validate_ipi_cell(payload.cell, + payload.inv_cell, + kMaxCellCondition, + kInverseAbsoluteTolerance, + kInverseRelativeTolerance); + 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(payload.nat_socket, + ucell.nat, + coordinate_count, + io_message)) + { + io_failed = 1; + } + if (io_failed == 0) + { + payload.positions = socket.read_doubles(coordinate_count); + if (!SocketFrame::validate_positions(payload.positions, + coordinate_count, + io_message)) + { + io_failed = 1; + } + } + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + } + quit_if_root_failed(io_failed, io_message); + return payload; +} + +void bcast_posdata(PosdataPayload& payload) +{ + bcast_int32(payload.nat_socket); + std::vector cell_values(payload.cell.begin(), payload.cell.end()); + std::vector inverse_values(payload.inv_cell.begin(), payload.inv_cell.end()); + bcast_double_vector(cell_values); + bcast_double_vector(inverse_values); + if (!is_root()) + { + payload.cell = {{0.0}}; + payload.inv_cell = {{0.0}}; + std::copy(cell_values.begin(), cell_values.end(), payload.cell.begin()); + std::copy(inverse_values.begin(), inverse_values.end(), payload.inv_cell.begin()); + if (payload.nat_socket >= 0) + { + payload.positions.assign(static_cast(3 * payload.nat_socket), 0.0); + } + } + bcast_double_vector(payload.positions); +} + +void check_posdata_geometry(const DriverContext& context, + const PosdataPayload& payload, + DriverContext& mutable_context) +{ + const double max_cell_delta_bohr = max_abs_delta( + std::vector(payload.cell.begin(), payload.cell.end()), + context.reference_cell); + if (max_cell_delta_bohr > unchanged_cell_tolerance(payload.cell)) + { + ModuleBase::WARNING_QUIT("ABACUS socket", + "variable-cell socket updates are not supported yet."); + } + if (!mutable_context.checked_initial_positions) + { + mutable_context.checked_initial_positions = true; + if (max_wrapped_direct_delta_from_unitcell(*context.ucell, payload.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."); + } + } +} + +void run_esolver_for_positions(UnitCell& ucell, + ModuleESolver::ESolver* esolver, + const std::vector& positions, + const int istep) +{ + 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 + { + 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"); + } +} + +void compute_energy_hartree(ModuleESolver::ESolver* esolver, + ComputedFrame& computed, + std::ofstream& ofs_running) +{ + double energy_ry = 0.0; + try + { + energy_ry = 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"); + } + const 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 * kRyToHartree; + if (is_root()) + { + ofs_running << " ABACUS socket return energy " + << energy_ry << " Ry, " + << energy_ry * ModuleBase::Ry_to_eV << " eV, " + << computed.energy_hartree << " Ha" << std::endl; + } +} + +void compute_forces(UnitCell& ucell, + ModuleESolver::ESolver* esolver, + ComputedFrame& computed) +{ + ModuleBase::matrix force; + try + { + 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"); + } + int local_failed = 0; + std::string local_message; + try + { + computed.forces_hartree_per_bohr = flatten_forces_hartree_per_bohr(force, ucell.nat); + } + 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; +} + +void compute_stress(UnitCell& ucell, + ModuleESolver::ESolver* esolver, + ComputedFrame& computed) +{ + ModuleBase::matrix stress; + try + { + 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"); + } + int local_failed = 0; + std::string local_message; + try + { + const SocketFrame::VirialConversion virial + = SocketFrame::make_ipi_virial(matrix9_from_stress(stress), + ucell.omega, + kStressAbsoluteTolerance, + kStressRelativeTolerance); + 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; +} +} // namespace + +void handle_posdata(IpiSocket& socket, + DriverContext& context, + std::ofstream& ofs_running) +{ + PosdataPayload payload = read_posdata(socket, *context.ucell, context.state); + bcast_posdata(payload); + check_posdata_geometry(context, payload, context); + run_esolver_for_positions(*context.ucell, context.esolver, payload.positions, context.istep); + + ComputedFrame computed; + computed.scf_converged = all_ranks_converged(context.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."); + } + compute_energy_hartree(context.esolver, computed, ofs_running); + if (context.inp->cal_force) + { + compute_forces(*context.ucell, context.esolver, computed); + } + if (context.inp->cal_stress) + { + compute_stress(*context.ucell, context.esolver, computed); + } + computed.valid = true; + context.published = computed; + ++context.istep; + context.state = DriverState::HasData; +} + +void handle_getforce(IpiSocket& socket, DriverContext& context) +{ + int io_failed = 0; + std::string io_message; + if (is_root()) + { + try + { + if (context.state != DriverState::HasData || !context.published.valid) + { + throw std::runtime_error("GETFORCE requires HAVEDATA state and a valid frame"); + } + socket.write_header("FORCEREADY"); + socket.write_double(context.published.energy_hartree); + socket.write_int32(static_cast(context.nat_return)); + const std::vector forces + = context.published.forces_present + ? context.published.forces_hartree_per_bohr + : std::vector(static_cast(3 * context.nat_return), 0.0); + socket.write_doubles(forces); + socket.write_doubles(vector_from_matrix9(context.published.virial_wire_hartree)); + const std::string extra = properties_extra(context.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) + { + io_failed = 1; + io_message = exc.what(); + } + } + quit_if_root_failed(io_failed, io_message); + context.published = ComputedFrame(); + context.state = DriverState::Ready; +} + +void handle_exit(std::ofstream& ofs_running) +{ + if (is_root()) + { + ofs_running << " ABACUS socket driver received i-PI EXIT" << std::endl; + } +} +} // namespace SocketHandlers diff --git a/source/source_relax/socket_handlers.h b/source/source_relax/socket_handlers.h new file mode 100644 index 00000000000..4681491733e --- /dev/null +++ b/source/source_relax/socket_handlers.h @@ -0,0 +1,54 @@ +#ifndef ABACUS_SOURCE_RELAX_SOCKET_HANDLERS_H +#define ABACUS_SOURCE_RELAX_SOCKET_HANDLERS_H + +#include "source_relax/socket_utils.h" +#include "source_relax/socket_ipi.h" + +#include +#include +#include + +class UnitCell; + +namespace ModuleESolver +{ +class ESolver; +} +class Input_para; + +namespace SocketHandlers +{ +struct DriverContext +{ + ModuleESolver::ESolver* esolver = nullptr; + UnitCell* ucell = nullptr; + const Input_para* inp = nullptr; + SocketUtils::DriverState state = SocketUtils::DriverState::NeedInit; + int istep = 0; + int nat_return = 0; + SocketUtils::ComputedFrame published; + std::vector reference_cell; + bool checked_initial_positions = false; +}; + +// Reads the next i-PI header on the root rank and broadcasts it to all ranks. +// Returns an empty string when the peer closed the connection while no frame +// is pending. Calls WARNING_QUIT on unrecoverable I/O failure. +std::string read_header_bcast(IpiSocket& socket, const SocketUtils::DriverState state); + +void handle_status(IpiSocket& socket, const SocketUtils::DriverState state); + +void handle_init(IpiSocket& socket, + DriverContext& context, + std::ofstream& ofs_running); + +void handle_posdata(IpiSocket& socket, + DriverContext& context, + std::ofstream& ofs_running); + +void handle_getforce(IpiSocket& socket, DriverContext& context); + +void handle_exit(std::ofstream& ofs_running); +} // namespace SocketHandlers + +#endif diff --git a/source/source_relax/socket_ipi.cpp b/source/source_relax/socket_ipi.cpp index d0ba19869aa..e84d8637ba2 100644 --- a/source/source_relax/socket_ipi.cpp +++ b/source/source_relax/socket_ipi.cpp @@ -18,7 +18,7 @@ static_assert(std::numeric_limits::is_iec559, namespace { -constexpr std::size_t IPI_HEADER_LEN = 12; +constexpr std::size_t kIpiHeaderLen = 12; std::string errno_message(const std::string& prefix) { @@ -27,7 +27,7 @@ std::string errno_message(const std::string& prefix) std::string trim_header(const char* data) { - std::string value(data, IPI_HEADER_LEN); + std::string value(data, kIpiHeaderLen); while (!value.empty() && value.back() == ' ') { value.pop_back(); @@ -37,12 +37,12 @@ std::string trim_header(const char* data) std::string padded_header(const std::string& header) { - if (header.size() > IPI_HEADER_LEN) + if (header.size() > kIpiHeaderLen) { throw std::runtime_error("i-PI header is longer than 12 bytes: " + header); } std::string out = header; - out.resize(IPI_HEADER_LEN, ' '); + out.resize(kIpiHeaderLen, ' '); return out; } @@ -145,7 +145,7 @@ void IpiSocket::close() std::string IpiSocket::read_header() { - char header[IPI_HEADER_LEN]; + char header[kIpiHeaderLen]; std::size_t done = 0; while (done < sizeof(header)) { diff --git a/source/source_relax/socket_utils.cpp b/source/source_relax/socket_utils.cpp new file mode 100644 index 00000000000..ccdfa38b369 --- /dev/null +++ b/source/source_relax/socket_utils.cpp @@ -0,0 +1,348 @@ +#include "socket_utils.h" + +#include "source_base/global_function.h" +#include "source_base/mathzone.h" +#include "source_base/parallel_common.h" +#include "source_base/parallel_reduce.h" +#include "source_base/timer.h" +#include "source_cell/unitcell.h" +#include "source_cell/update_cell.h" + +#ifdef __MPI +#include +#endif + +#include +#include +#include +#include +#include +#include +#include + +namespace SocketUtils +{ +bool all_ranks_converged(const bool local_converged) +{ + int converged = local_converged ? 1 : 0; + Parallel_Reduce::reduce_min(converged); + return converged != 0; +} + +void throw_if_any_rank_failed(int local_failed, std::string local_message) +{ + int any_failed = local_failed; + Parallel_Reduce::reduce_max(any_failed); + if (any_failed != 0) + { + if (local_message.empty()) + { + local_message = "socket frame validation failed on another MPI rank"; + } + throw std::runtime_error(local_message); + } +} + +/// @brief Terminate the socket run after a failure inside a collective stage. +/// +/// Under MPI the failure may deadlock other ranks waiting in a collective, +/// so print a diagnostic and abort the whole job. In serial builds there is +/// no collective: the exception propagates to Socket_Driver::socket_driver, +/// which reports a clean WARNING_QUIT into the log. Do not replace the +/// serial throw with abort(); serial builds are user-facing. +[[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() +{ + int rank = kIpiRankRoot; +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &rank); +#endif + return rank == kIpiRankRoot; +} + +void bcast_double_vector(std::vector& values) +{ + if (!values.empty()) + { + Parallel_Common::bcast_double(values.data(), static_cast(values.size())); + } +} + +void bcast_int(int& value) +{ + Parallel_Common::bcast_int(value); +} + +void bcast_int32(std::int32_t& value) +{ + int tmp = static_cast(value); + Parallel_Common::bcast_int(tmp); + value = static_cast(tmp); +} + +void bcast_chars(char* value, const int size) +{ + Parallel_Common::bcast_char(value, size); +} + +void bcast_string(std::string& value) +{ + int size = static_cast(value.size()); + bcast_int(size); + if (!is_root()) + { + value.resize(static_cast(size)); + } + if (size > 0) + { + bcast_chars(&value[0], size); + } +} + +void quit_if_root_failed(int root_failed, std::string root_message) +{ + bcast_int(root_failed); + bcast_string(root_message); + if (root_failed != 0) + { + ModuleBase::WARNING_QUIT("ABACUS socket", root_message.empty() ? "i-PI socket I/O failed" : root_message); + } +} + +std::string bcast_header(std::string header) +{ + bcast_string(header); + return header; +} + +std::string address() +{ + const char* env = std::getenv("ABACUS_SOCKET_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_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()) + { + 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; +} + +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)) + { + ModuleBase::WARNING_QUIT("ABACUS socket", "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, const int nat) +{ + if (nat < 0 || force.nr != nat || 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) + { + 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 * kRyToHartree; + } + } + 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 SocketUtils diff --git a/source/source_relax/socket_utils.h b/source/source_relax/socket_utils.h new file mode 100644 index 00000000000..ccac9761967 --- /dev/null +++ b/source/source_relax/socket_utils.h @@ -0,0 +1,71 @@ +#ifndef ABACUS_SOURCE_RELAX_SOCKET_UTILS_H +#define ABACUS_SOURCE_RELAX_SOCKET_UTILS_H + +#include "source_relax/socket_frame.h" +#include "source_base/matrix.h" + +#include +#include +#include +#include + +class UnitCell; + +namespace SocketUtils +{ +constexpr double kRyToHartree = 0.5; +constexpr int kIpiRankRoot = 0; +constexpr double kMaxCellCondition = 1.0e12; +constexpr double kInverseAbsoluteTolerance + = 64.0 * std::numeric_limits::epsilon(); +constexpr double kInverseRelativeTolerance = 64.0; +constexpr double kStressAbsoluteTolerance = 1.0e-10; +constexpr double kStressRelativeTolerance = 1.0e-8; +constexpr std::int32_t kMaxInitBytes = static_cast(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(bool local_converged); +void throw_if_any_rank_failed(int local_failed, std::string local_message); +[[noreturn]] void fail_during_collective_stage(const char* stage, + const std::string& message); +std::string properties_extra(const ComputedFrame& frame); +bool is_root(); +void bcast_double_vector(std::vector& values); +void bcast_int(int& value); +void bcast_int32(std::int32_t& value); +void bcast_chars(char* value, int size); +void bcast_string(std::string& value); +void quit_if_root_failed(int root_failed, std::string root_message); +std::string bcast_header(std::string header); +std::string address(); +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); +double max_abs_delta(const std::vector& a, const std::vector& b); +double unchanged_cell_tolerance(const SocketFrame::Matrix9& cell); +void set_positions_from_ipi_bohr(UnitCell& ucell, + const std::vector& positions_bohr); +std::vector flatten_forces_hartree_per_bohr(const ModuleBase::matrix& force, + int nat); +SocketFrame::Matrix9 matrix9_from_stress(const ModuleBase::matrix& stress); +std::vector vector_from_matrix9(const SocketFrame::Matrix9& values); +} // namespace SocketUtils + +#endif diff --git a/source/source_relax/test/CMakeLists.txt b/source/source_relax/test/CMakeLists.txt index a5fdadb78c1..22376932aa3 100644 --- a/source/source_relax/test/CMakeLists.txt +++ b/source/source_relax/test/CMakeLists.txt @@ -9,19 +9,23 @@ install(DIRECTORY support DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) AddTest( TARGET MODULE_RELAX_socket_ipi_test - SOURCES socket_ipi_test.cpp ../socket_ipi.cpp + SOURCES test_socket_ipi.cpp ../socket_ipi.cpp ) AddTest( TARGET MODULE_RELAX_socket_frame_test - SOURCES socket_frame_test.cpp ../socket_frame.cpp + SOURCES test_socket_frame.cpp ../socket_frame.cpp ../socket_frame_utils.cpp ) AddTest( TARGET MODULE_RELAX_socket_driver_test - LIBS base device - SOURCES socket_driver_test.cpp + LIBS MPI::MPI_CXX base device + KEEP_FEATURE_DEFINITIONS __MPI + SOURCES test_socket_driver.cpp ../socket_driver.cpp + ../socket_utils.cpp + ../socket_handlers.cpp + ../socket_frame_utils.cpp ../socket_frame.cpp ../socket_ipi.cpp ../../source_cell/update_cell.cpp diff --git a/source/source_relax/test/for_test.h b/source/source_relax/test/for_test.h index ed16ca04cc4..043c88ed470 100644 --- a/source/source_relax/test/for_test.h +++ b/source/source_relax/test/for_test.h @@ -84,6 +84,14 @@ Atom::Atom() Atom::~Atom() { } +#ifdef __MPI +void Atom::bcast_atom() +{ +} +void Atom::bcast_atom2() +{ +} +#endif Atom_pseudo::Atom_pseudo() { } diff --git a/source/source_relax/test/socket_driver_test.cpp b/source/source_relax/test/test_socket_driver.cpp similarity index 98% rename from source/source_relax/test/socket_driver_test.cpp rename to source/source_relax/test/test_socket_driver.cpp index fe361e55d65..5501ef9475d 100644 --- a/source/source_relax/test/socket_driver_test.cpp +++ b/source/source_relax/test/test_socket_driver.cpp @@ -2,6 +2,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include "mpi.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" #include "source_io/module_parameter/input_parameter.h" @@ -26,7 +27,7 @@ namespace { -constexpr std::size_t IPI_HEADER_LEN = 12; +constexpr std::size_t kIpiHeaderLen = 12; std::string errno_message(const std::string& prefix) { @@ -70,7 +71,7 @@ void send_value(const int fd, const T& value) void send_header(const int fd, const std::string& header) { std::string padded = header; - padded.resize(IPI_HEADER_LEN, ' '); + padded.resize(kIpiHeaderLen, ' '); send_all(fd, padded.data(), padded.size()); } @@ -93,7 +94,7 @@ bool try_send_status(const int fd) std::string read_header_or_close(const int fd) { - char header[IPI_HEADER_LEN]; + char header[kIpiHeaderLen]; std::size_t done = 0; while (done < sizeof(header)) { @@ -613,3 +614,12 @@ TEST(SocketDriverTest, InvalidNextFrameCannotReturnPreviousResults) EXPECT_NE(0, result.exit_code); EXPECT_THAT(result.diagnostic, testing::HasSubstr("finite")); } + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); + MPI_Finalize(); + return result; +} diff --git a/source/source_relax/test/socket_frame_test.cpp b/source/source_relax/test/test_socket_frame.cpp similarity index 100% rename from source/source_relax/test/socket_frame_test.cpp rename to source/source_relax/test/test_socket_frame.cpp diff --git a/source/source_relax/test/socket_ipi_test.cpp b/source/source_relax/test/test_socket_ipi.cpp similarity index 98% rename from source/source_relax/test/socket_ipi_test.cpp rename to source/source_relax/test/test_socket_ipi.cpp index c2b7f415334..2f02e1cd55a 100644 --- a/source/source_relax/test/socket_ipi_test.cpp +++ b/source/source_relax/test/test_socket_ipi.cpp @@ -18,7 +18,7 @@ namespace { -constexpr std::size_t IPI_HEADER_LEN = 12; +constexpr std::size_t kIpiHeaderLen = 12; std::string errno_message(const std::string& prefix) { @@ -74,7 +74,7 @@ void recv_all(int fd, void* data, std::size_t nbytes) std::string padded_header(const std::string& header) { std::string padded = header; - padded.resize(IPI_HEADER_LEN, ' '); + padded.resize(kIpiHeaderLen, ' '); return padded; } @@ -170,7 +170,7 @@ TEST(IpiSocketTest, WriteHeaderPadsToTwelveBytes) try { const int fd = server.accept_once(); - char buffer[IPI_HEADER_LEN]; + char buffer[kIpiHeaderLen]; recv_all(fd, buffer, sizeof(buffer)); received.assign(buffer, sizeof(buffer)); ::close(fd);