diff --git a/python/pyabacus/src/ModuleDriver/py_driver.cpp b/python/pyabacus/src/ModuleDriver/py_driver.cpp index 26b02bf49f8..2083b997127 100644 --- a/python/pyabacus/src/ModuleDriver/py_driver.cpp +++ b/python/pyabacus/src/ModuleDriver/py_driver.cpp @@ -451,7 +451,7 @@ CalculationResult PyDriver::run(const std::string& input_dir, if (cal == "scf" || cal == "relax" || cal == "cell-relax" || cal == "nscf") { Relax_Driver rl_driver; - rl_driver.relax_driver(impl_->p_esolver_, *impl_->ucell_, PARAM.inp, GlobalV::ofs_running); + rl_driver.relax_driver(impl_->p_esolver_, *impl_->ucell_, PARAM, GlobalV::ofs_running); } else if (cal == "get_s") { diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 142750d9fb3..f5705420ba2 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -271,6 +271,7 @@ OBJS_DEEPKS=lcao_deepks.o\ OBJS_ELECSTAT=elecstate.o\ + electronic_occupations.o\ init_scf.o\ param_update.o\ estate_e_terms.o\ @@ -298,6 +299,7 @@ OBJS_ELECSTAT=elecstate.o\ efield.o\ gatefield.o\ potential_new.o\ + electrostatic_potential.o\ potential_types.o\ pot_surchem.o\ pot_sep.o\ @@ -627,6 +629,7 @@ OBJS_RELAXATION=relax_data.o\ socket_driver.o\ socket_utils.o\ socket_handlers.o\ + socket_constant_potential.o\ cg_base.o\ bfgs_basic.o\ relax_driver.o\ diff --git a/source/source_esolver/esolver.h b/source/source_esolver/esolver.h index f302584cbb8..b89ec6bac66 100644 --- a/source/source_esolver/esolver.h +++ b/source/source_esolver/esolver.h @@ -6,6 +6,14 @@ #include "source_cell/unitcell.h" struct Input_para; +class Parameter; + +namespace elecstate +{ +struct ElectronicControl; +struct PotentialWindow; +struct ElectronicObservables; +} namespace ModuleESolver { @@ -46,6 +54,11 @@ class ESolver //! calcualte stress of given cell virtual void cal_stress(BaseCell& cell, ModuleBase::matrix& stress) = 0; + // Optional first-principles capability. Unsupported solvers return false. + virtual bool apply_electronic_control(Parameter&, const elecstate::ElectronicControl&) { return false; } + virtual bool electronic_observables(const UnitCell&, const Parameter&, const elecstate::PotentialWindow&, + elecstate::ElectronicObservables&) const { return false; } + bool conv_esolver = true; // whether esolver is converged std::string classname; diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 070e860f388..554d69b847f 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -1,4 +1,10 @@ #include "esolver_fp.h" +#include "source_estate/electronic_control.h" +#include "source_estate/electronic_occupations.h" +#include "source_estate/elecstate_tools.h" +#include "source_estate/module_pot/electrostatic_potential.h" + +#include #include "source_basis/module_pw/pw_basis_big.h" // use PW_Basis_Big #include "source_base/tool_quit.h" @@ -24,6 +30,34 @@ #include "source_base/module_parallel/para_bridge.h" namespace ModuleESolver { +bool ESolver_FP::apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) +{ + if (!pelec || !pelec->charge || !pelec->pot) + { throw std::runtime_error("electronic state is not initialized for electron control"); } + elecstate::ParamUpdater::update_electron_count(control.nelect, parameters); + elecstate::update_electron_count(*pelec, control.nelect, parameters.inp.nspin, + parameters.inp.nupdown, parameters.globalv.two_fermi); + return true; +} + +bool ESolver_FP::electronic_observables(const UnitCell& ucell, const Parameter& parameters, + const elecstate::PotentialWindow& window, + elecstate::ElectronicObservables& result) const +{ + if (!this->conv_esolver || !pelec || !pelec->charge || !pelec->pot) + { throw std::runtime_error("electronic observables require a converged SCF state"); } + const int npool = parameters.inp.kpar * parameters.inp.bndpar; + result.nelect = elecstate::checked_occupation_electron_count(pelec->wg, parameters.inp.nelec, + npool, pw_rhod->poolnproc); + result.fermi_ry = elecstate::electronic_chemical_potential(pelec->eferm, parameters.globalv.two_fermi); + if (window.axis >= 0) + { + result.reference_ry = elecstate::electrostatic_reference(ucell, *pw_rhod, *pelec->charge, + *pelec->pot, solvent, parameters.inp, window); + } + return true; +} + ESolver_FP::ESolver_FP() { diff --git a/source/source_esolver/esolver_fp.h b/source/source_esolver/esolver_fp.h index 39612979a9a..f9c39d3788b 100644 --- a/source/source_esolver/esolver_fp.h +++ b/source/source_esolver/esolver_fp.h @@ -41,6 +41,10 @@ class ESolver_FP : public ESolver virtual void after_all_runners(BaseCell& basecell) override; + bool apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) override; + bool electronic_observables(const UnitCell& ucell, const Parameter& parameters, + const elecstate::PotentialWindow& window, elecstate::ElectronicObservables& result) const override; + protected: virtual void before_scf(UnitCell& ucell, const int istep); diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index dcb2a3c61f6..26170a09ad9 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -179,12 +179,14 @@ void ESolver_KS::runner(BaseCell& basecell, const int istep) ModuleBase::TITLE("ESolver_KS", "runner"); ModuleBase::timer::start(this->classname, "runner"); + this->conv_esolver = false; + // 1) before_scf (electronic iteration loops) this->before_scf(ucell, istep); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF"); const int iter_limit = this->scf_iter_limit(istep); - const bool converged = this->scf_loop(ucell, istep, iter_limit); + const bool converged = this->scf_loop(ucell, istep, iter_limit) && !this->oscillate_esolver; this->after_scf(ucell, istep, converged); this->conv_esolver = converged; ModuleBase::timer::end(this->classname, "runner"); diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index f476cf78d90..2e0b1c7faa5 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -162,6 +162,17 @@ void ESolver_KS_LCAO::before_all_runners(BaseCell& basecell, const Input return; } +template +bool ESolver_KS_LCAO::apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) +{ + if (!ESolver_FP::apply_electronic_control(parameters, control)) { return false; } +#ifdef __EXX + if (this->exx_nao.exd) { this->exx_nao.exd->begin_electronic_step(); } + if (this->exx_nao.exc) { this->exx_nao.exc->begin_electronic_step(); } +#endif + return true; +} + template void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { diff --git a/source/source_esolver/esolver_ks_lcao.h b/source/source_esolver/esolver_ks_lcao.h index 6fca9c417ec..032facf4207 100644 --- a/source/source_esolver/esolver_ks_lcao.h +++ b/source/source_esolver/esolver_ks_lcao.h @@ -35,6 +35,8 @@ class ESolver_KS_LCAO : public ESolver_KS //! Release the electronic wave functions. ~ESolver_KS_LCAO(); + bool apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) override; + //! Initialize the LCAO basis, electronic state, and potentials before all steps. void before_all_runners(BaseCell& basecell, const Input_para& inp) override; diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index c700eefa055..330e0193784 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -18,6 +18,7 @@ #include "source_estate/module_charge/chg_routine.h" // use charge mixing, mohan add 20251006 #include "source_estate/setup_estate_pw.h" // mohan add 20251005 +#include "source_hamilt/module_xc/xc_functional.h" #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type used via general_exx_info_ #include "source_io/module_ctrl/ctrl_output_pw.h" // mohan add 20250927 #include "source_pw/module_pwdft/dftu_base.h" @@ -159,6 +160,16 @@ void ESolver_KS_PW::before_all_runners(BaseCell& basecell, const Inpu this->exx_helper->init(ucell, inp, this->pelec->wg, this->general_exx_info_); } +template +bool ESolver_KS_PW::apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) +{ + const int functional_type = XC_Functional::get_func_type(); + validate_electronic_step_exx(functional_type, this->general_exx_info_.cal_exx); + if (!ESolver_FP::apply_electronic_control(parameters, control)) { return false; } + this->exx_helper->begin_electronic_step(); + return true; +} + template void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) { diff --git a/source/source_esolver/esolver_ks_pw.h b/source/source_esolver/esolver_ks_pw.h index a03e73c943e..53257a91c36 100644 --- a/source/source_esolver/esolver_ks_pw.h +++ b/source/source_esolver/esolver_ks_pw.h @@ -34,6 +34,8 @@ class ESolver_KS_PW : public ESolver_KS Exx_HelperBase* exx_helper = nullptr; + bool apply_electronic_control(Parameter& parameters, const elecstate::ElectronicControl& control) override; + protected: virtual void before_scf(UnitCell& ucell, const int istep) override; diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index 9f83433a275..5ff56c6bba8 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -2,6 +2,7 @@ list( APPEND objects elecstate.cpp + electronic_occupations.cpp init_scf.cpp estate_e_terms.cpp elecstate_energy.cpp @@ -19,6 +20,7 @@ list( module_pot/gatefield.cpp module_pot/efield.cpp module_pot/h_hartree_pw.cpp + module_pot/electrostatic_potential.cpp module_pot/h_tddft_pw.cpp module_pot/pot_xc.cpp module_pot/pot_local.cpp diff --git a/source/source_estate/elecstate.cpp b/source/source_estate/elecstate.cpp index 5ed97e678c7..7b012a3c8f7 100644 --- a/source/source_estate/elecstate.cpp +++ b/source/source_estate/elecstate.cpp @@ -16,11 +16,16 @@ const double* ElecState::getRho(int spin) const void ElecState::init_nelec_spin() { - this->nelec_spin.resize(PARAM.inp.nspin); - if (PARAM.inp.nspin == 2) + this->init_nelec_spin(PARAM.inp.nspin, PARAM.inp.nelec, PARAM.inp.nupdown); +} + +void ElecState::init_nelec_spin(int nspin, double nelect, double nupdown) +{ + this->nelec_spin.resize(nspin); + if (nspin == 2) { - this->nelec_spin[0] = (PARAM.inp.nelec + PARAM.inp.nupdown) / 2.0; - this->nelec_spin[1] = (PARAM.inp.nelec - PARAM.inp.nupdown) / 2.0; + this->nelec_spin[0] = (nelect + nupdown) / 2.0; + this->nelec_spin[1] = (nelect - nupdown) / 2.0; } } diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index ba96ff14576..a5643e6d7cc 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -116,6 +116,7 @@ class ElecState // if nupdown is not 0(TWO_EFERMI case), // nelec_spin will be fixed and weights will be constrained void init_nelec_spin(); + void init_nelec_spin(int nspin, double nelect, double nupdown); // used to record number of electrons per spin index // for NSPIN=2, it will record number of spin up and number of spin down // for NSPIN=4, it will record total number, magnetization for x, y, z direction diff --git a/source/source_estate/elecstate_tools.cpp b/source/source_estate/elecstate_tools.cpp index 3b64aca047f..c1834c18fa9 100644 --- a/source/source_estate/elecstate_tools.cpp +++ b/source/source_estate/elecstate_tools.cpp @@ -1,6 +1,9 @@ #include "elecstate_tools.h" #include "occupy.h" + +#include +#include #include "source_base/module_parallel/para_band_output.h" #include "source_base/module_parallel/para_bridge.h" #include "source_base/parallel_reduce.h" @@ -9,6 +12,25 @@ namespace elecstate { +void update_electron_count(ElecState& state, double nelect, int nspin, double nupdown, bool two_fermi) +{ + state.eferm.two_efermi = two_fermi; + state.init_nelec_spin(nspin, nelect, nupdown); +} + +double electronic_chemical_potential(const Efermi& fermi, bool two_fermi) +{ + double chemical_potential = fermi.ef; + if (two_fermi) + { + // Total-N derivative at fixed N_up-N_down. + chemical_potential = 0.5 * (fermi.ef_up + fermi.ef_dw); + } + if (!std::isfinite(chemical_potential)) + { throw std::runtime_error("non-finite SCF chemical potential"); } + return chemical_potential; +} + void calEBand(const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg, fenergy& f_en) { ModuleBase::TITLE("ElecState", "calEBand"); diff --git a/source/source_estate/elecstate_tools.h b/source/source_estate/elecstate_tools.h index bf1f4bb4274..d7471cb4d80 100644 --- a/source/source_estate/elecstate_tools.h +++ b/source/source_estate/elecstate_tools.h @@ -5,6 +5,9 @@ namespace elecstate { +void update_electron_count(ElecState& state, double nelect, int nspin, double nupdown, bool two_fermi); +double electronic_chemical_potential(const Efermi& fermi, bool two_fermi); + void calEBand(const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg, fenergy& f_en); void calculate_weights(const ModuleBase::matrix& ekb, diff --git a/source/source_estate/electronic_control.h b/source/source_estate/electronic_control.h new file mode 100644 index 00000000000..9bd7839fe97 --- /dev/null +++ b/source/source_estate/electronic_control.h @@ -0,0 +1,27 @@ +#ifndef ABACUS_ELECTRONIC_CONTROL_H +#define ABACUS_ELECTRONIC_CONTROL_H + +namespace elecstate +{ +// Physical requests and results; no socket or serialization dependencies. +struct ElectronicControl +{ + double nelect = 0.0; +}; + +struct PotentialWindow +{ + int axis = -1; + double minimum_A = 0.0; + double maximum_A = 0.0; +}; + +struct ElectronicObservables +{ + double nelect = 0.0; + double fermi_ry = 0.0; + double reference_ry = 0.0; +}; +} + +#endif diff --git a/source/source_estate/electronic_occupations.cpp b/source/source_estate/electronic_occupations.cpp new file mode 100644 index 00000000000..37917edabf1 --- /dev/null +++ b/source/source_estate/electronic_occupations.cpp @@ -0,0 +1,65 @@ +#include "electronic_occupations.h" +#include "source_base/matrix.h" +#include "source_base/parallel_reduce.h" + +#include +#include +#include +#include + +namespace elecstate +{ +double checked_occupation_electron_count(const ModuleBase::matrix& weights, + double requested_nelect, + int npool, + int ranks_per_pool) +{ + int invalid = !std::isfinite(requested_nelect) || requested_nelect <= 0.0 + || npool <= 0 || ranks_per_pool <= 0; + double occupied = 0.0; + double compensation = 0.0; + for (int ik = 0; ik < weights.nr; ++ik) + { + for (int ib = 0; ib < weights.nc; ++ib) + { + const double weight = weights(ik, ib); + if (!std::isfinite(weight)) + { + invalid = 1; + continue; + } + // Compensated sum of the local k/spin and band shard. Generalized + // smearing weights may be negative; only their total is constrained. + const double corrected = weight - compensation; + const double next = occupied + corrected; + compensation = (next - occupied) - corrected; + occupied = next; + } + } + if (!std::isfinite(occupied)) { invalid = 1; } + // Synchronize failures before a rank can enter a different collective path. + Parallel_Reduce::reduce_all(invalid); + if (invalid) + { + throw std::runtime_error("constant-potential electron count requires finite occupations " + "and a valid target/occupation partition"); + } + + // Same distribution as the native band-energy sum: remove the replicated + // copies inside each pool before combining KPAR and BNDPAR contributions. + Parallel_Reduce::reduce_double_allpool(npool, ranks_per_pool, occupied); + const double tolerance = 1.0e-8; + int mismatch = !std::isfinite(occupied) || occupied <= 0.0 + || std::abs(occupied - requested_nelect) > tolerance; + Parallel_Reduce::reduce_all(mismatch); + if (mismatch) + { + std::ostringstream message; + message << std::setprecision(17) + << "constant-potential occupation electron count mismatch: requested=" << requested_nelect + << ", occupied=" << occupied << ", absolute tolerance=" << tolerance; + throw std::runtime_error(message.str()); + } + return occupied; +} +} diff --git a/source/source_estate/electronic_occupations.h b/source/source_estate/electronic_occupations.h new file mode 100644 index 00000000000..3b014171b60 --- /dev/null +++ b/source/source_estate/electronic_occupations.h @@ -0,0 +1,20 @@ +#ifndef ABACUS_ELECTRONIC_OCCUPATIONS_H +#define ABACUS_ELECTRONIC_OCCUPATIONS_H + +namespace ModuleBase +{ +class matrix; +} + +namespace elecstate +{ +// weights already include k-point weights and spin degeneracy. Rows/bands are +// partitioned across npool groups and replicated on ranks_per_pool ranks. +// All ranks must participate; invalid occupations or a mismatch fail together. +double checked_occupation_electron_count(const ModuleBase::matrix& weights, + double requested_nelect, + int npool, + int ranks_per_pool); +} + +#endif diff --git a/source/source_estate/module_pot/electrostatic_potential.cpp b/source/source_estate/module_pot/electrostatic_potential.cpp new file mode 100644 index 00000000000..933ebb7ee12 --- /dev/null +++ b/source/source_estate/module_pot/electrostatic_potential.cpp @@ -0,0 +1,106 @@ +#include "electrostatic_potential.h" +#include "potential_new.h" +#include "source_estate/electronic_control.h" +#include "h_hartree_pw.h" +#include "efield.h" +#include "source_base/module_parallel/para_mpi_func.h" +#include "source_base/module_parallel/para_world.h" +#include "source_estate/module_charge/charge.h" +#include "source_hamilt/module_surchem/surchem.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include + +namespace elecstate +{ +double electrostatic_reference(const UnitCell& cell, ModulePW::PW_Basis& basis, + const Charge& charge, const Potential& pot, + const surchem& solvent, const Input_para& input, + const PotentialWindow& window) +{ + const double* fixed = basis.nrxx > 0 ? pot.get_fixed_v() : nullptr; + const auto* pcc_potential = pot.pcc_electron_potential(); + const auto* solvent_potential = pot.find_sol_pot(); + const auto potential = electrostatic_potential(cell, basis, charge, fixed, solvent, input, + pcc_potential, solvent_potential); + return planar_potential_average(cell, basis, potential, window.axis, window.minimum_A, window.maximum_A); +} + +std::vector electrostatic_potential(const UnitCell& cell, ModulePW::PW_Basis& basis, + const Charge& charge, const double* fixed, + const surchem& solvent, const Input_para& input, + const std::vector* pcc_potential, + const std::vector* solvent_potential) +{ + if (basis.nrxx > 0 && (!fixed || (input.imp_sol == 1 && !solvent.delta_phi))) + { throw std::runtime_error("electrostatic reference requires initialized fixed and solvent potentials"); } + if (solvent_potential && solvent_potential->size() != static_cast(basis.nrxx)) + { throw std::runtime_error("Solvent potential does not match the output grid"); } + const auto hartree = H_Hartree_pw::v_hartree(cell, &basis, input.nspin, charge.rho); + ModuleBase::matrix dipole; + const bool add_dipole = input.efield_flag && input.dip_cor_flag; + if (add_dipole) + { dipole = Efield::add_efield(cell, &basis, input.nspin, charge.rho, solvent); } + std::vector potential(basis.nrxx); + for (int ir = 0; ir < basis.nrxx; ++ir) + { + potential[ir] = hartree(0, ir) + fixed[ir]; + if (add_dipole) { potential[ir] += dipole(0, ir); } + if (pcc_potential) { potential[ir] += (*pcc_potential)[ir]; } + if (solvent_potential) { potential[ir] += (*solvent_potential)[ir]; } + if (input.imp_sol == 1) { potential[ir] += solvent.delta_phi[ir]; } + } + return potential; +} + +double planar_potential_average(const UnitCell& cell, const ModulePW::PW_Basis& basis, + const std::vector& potential, int axis, + double minimum_A, double maximum_A) +{ + if (axis < 0 || axis > 2 || !std::isfinite(minimum_A) || !std::isfinite(maximum_A) + || minimum_A < 0 || maximum_A <= minimum_A) + { throw std::runtime_error("invalid electrostatic averaging axis/window"); } + const double lattice[3][3] = { + {cell.latvec.e11, cell.latvec.e12, cell.latvec.e13}, + {cell.latvec.e21, cell.latvec.e22, cell.latvec.e23}, + {cell.latvec.e31, cell.latvec.e32, cell.latvec.e33}}; + const int first = (axis + 1) % 3; + const int second = (axis + 2) % 3; + const double cross[3] = { + lattice[first][1] * lattice[second][2] - lattice[first][2] * lattice[second][1], + lattice[first][2] * lattice[second][0] - lattice[first][0] * lattice[second][2], + lattice[first][0] * lattice[second][1] - lattice[first][1] * lattice[second][0]}; + double area2 = 0.0; + double volume = 0.0; + for (int i = 0; i < 3; ++i) { area2 += cross[i]*cross[i]; volume += lattice[axis][i]*cross[i]; } + const double height_A = cell.lat0_angstrom * std::abs(volume) / std::sqrt(area2); + if (!std::isfinite(height_A) || height_A <= 0.0 || maximum_A > height_A + 1.0e-10) + { throw std::runtime_error("electrostatic averaging window lies outside the cell height"); } + const int grid[3] = {basis.nx, basis.ny, basis.nz}; + if (grid[0] <= 0 || grid[1] <= 0 || grid[2] <= 0 || basis.nrxx < 0 + || (basis.nrxx > 0 && basis.nplane <= 0) || potential.size() != static_cast(basis.nrxx)) + { throw std::runtime_error("invalid distributed FFT grid for electrostatic averaging"); } + double reduction[3] = {0.0, 0.0, 0.0}; + for (int ir = 0; ir < basis.nrxx; ++ir) + { + const int xy = ir / basis.nplane; + const int indices[3] = {xy / basis.ny, xy % basis.ny, ir % basis.nplane + basis.startz_current}; + // FFT points and ABACUS cube origins are at i/n, not voxel centers. + const double coordinate_A = double(indices[axis]) / grid[axis] * height_A; + if (!std::isfinite(potential[ir])) { reduction[2] += 1.0; } + if (coordinate_A >= minimum_A && coordinate_A < maximum_A) + { reduction[0] += potential[ir]; reduction[1] += 1.0; } + } + // Reduce over this FFT grid's communicator without summing replicated pools. +#ifdef __MPI + const auto grid_world = Parallel::ParaWorld::make_mpi("electrostatic_grid", basis.pool_world); +#else + const auto grid_world = Parallel::ParaWorld::serial("electrostatic_grid"); +#endif + Parallel::reduce_all(reduction, 3, grid_world); + if (reduction[2] || reduction[1] <= 0.0 || !std::isfinite(reduction[0])) + { throw std::runtime_error("non-finite electrostatic potential or no FFT points in averaging window"); } + return reduction[0] / reduction[1]; +} +} diff --git a/source/source_estate/module_pot/electrostatic_potential.h b/source/source_estate/module_pot/electrostatic_potential.h new file mode 100644 index 00000000000..64dcbc43bee --- /dev/null +++ b/source/source_estate/module_pot/electrostatic_potential.h @@ -0,0 +1,30 @@ +#ifndef ABACUS_ELECTROSTATIC_POTENTIAL_H +#define ABACUS_ELECTROSTATIC_POTENTIAL_H + +#include +class UnitCell; +class Charge; +class surchem; +struct Input_para; +namespace ModulePW { class PW_Basis; } + +namespace elecstate +{ +class Potential; +struct PotentialWindow; +double electrostatic_reference(const UnitCell& cell, ModulePW::PW_Basis& basis, + const Charge& charge, const Potential& pot, + const surchem& solvent, const Input_para& input, + const PotentialWindow& window); + +// Same local electrostatic reference for cube output and electronic observables. +std::vector electrostatic_potential(const UnitCell& cell, ModulePW::PW_Basis& basis, + const Charge& charge, const double* fixed, + const surchem& solvent, const Input_para& input, + const std::vector* pcc_potential, + const std::vector* solvent_potential); +double planar_potential_average(const UnitCell& cell, const ModulePW::PW_Basis& basis, + const std::vector& potential, int axis, + double minimum_A, double maximum_A); +} +#endif diff --git a/source/source_estate/param_update.cpp b/source/source_estate/param_update.cpp index 3ff80949651..90c5c30cc7a 100644 --- a/source/source_estate/param_update.cpp +++ b/source/source_estate/param_update.cpp @@ -14,4 +14,9 @@ void ParamUpdater::update_from_atoms_info(const AtomsInfoResult& atoms_info) PARAM.sys.ks_run = atoms_info.ks_run; } -} \ No newline at end of file +void ParamUpdater::update_electron_count(double nelect, Parameter& parameters) +{ + parameters.input.nelec = nelect; +} + +} diff --git a/source/source_estate/param_update.h b/source/source_estate/param_update.h index 89a0d90023b..d63f3870584 100644 --- a/source/source_estate/param_update.h +++ b/source/source_estate/param_update.h @@ -3,13 +3,16 @@ #include "source_cell/cal_atoms_info.h" +class Parameter; + namespace elecstate { class ParamUpdater { public: static void update_from_atoms_info(const AtomsInfoResult& atoms_info); + static void update_electron_count(double nelect, Parameter& parameters); }; } -#endif \ No newline at end of file +#endif diff --git a/source/source_estate/unittests/CMakeLists.txt b/source/source_estate/unittests/CMakeLists.txt index 7cc711ec56e..bdb67dff6a1 100644 --- a/source/source_estate/unittests/CMakeLists.txt +++ b/source/source_estate/unittests/CMakeLists.txt @@ -26,6 +26,41 @@ addtest( if(ENABLE_MPI) + addtest( + TARGET MODULE_ESTATE_electronic_occupations_test + LIBS base device MPI::MPI_CXX + KEEP_FEATURE_DEFINITIONS __MPI + SOURCES test_electronic_occupations.cpp ../electronic_occupations.cpp) + set_tests_properties(MODULE_ESTATE_electronic_occupations_test + PROPERTIES TIMEOUT 30 ENVIRONMENT "OMP_NUM_THREADS=1") + add_test( + NAME MODULE_ESTATE_electronic_occupations_mpi4_test + COMMAND + ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + $ + ${MPIEXEC_POSTFLAGS}) + set_tests_properties( + MODULE_ESTATE_electronic_occupations_mpi4_test + PROPERTIES TIMEOUT 30 PROCESSORS 4 ENVIRONMENT "OMP_NUM_THREADS=1") + + addtest( + TARGET MODULE_ESTATE_electrostatic_potential_test + LIBS parameter base device planewave MPI::MPI_CXX + KEEP_FEATURE_DEFINITIONS __MPI + SOURCES test_electrostatic_potential.cpp + ../module_pot/electrostatic_potential.cpp) + set_tests_properties(MODULE_ESTATE_electrostatic_potential_test + PROPERTIES TIMEOUT 30 ENVIRONMENT "OMP_NUM_THREADS=1") + add_test( + NAME MODULE_ESTATE_electrostatic_potential_mpi4_test + COMMAND + ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + $ + ${MPIEXEC_POSTFLAGS}) + set_tests_properties( + MODULE_ESTATE_electrostatic_potential_mpi4_test + PROPERTIES TIMEOUT 30 PROCESSORS 4 ENVIRONMENT "OMP_NUM_THREADS=1") + set(pcc_test_sources ../module_pot/sccs_output.cpp ../module_pot/pot_pcc.cpp diff --git a/source/source_estate/unittests/test_elecstate.cpp b/source/source_estate/unittests/test_elecstate.cpp index 317ef2805a2..5d157875248 100644 --- a/source/source_estate/unittests/test_elecstate.cpp +++ b/source/source_estate/unittests/test_elecstate.cpp @@ -213,6 +213,32 @@ TEST_F(ElecStateTest, InitNelecSpin) EXPECT_EQ(elecstate->nelec_spin[1], 5.0); } +TEST_F(ElecStateTest, ElectronicRequestUsesExplicitSpinState) +{ + // The fixture uses global N=10, nspin=1 and two_fermi=false. + elecstate::update_electron_count(*elecstate, 8.2, 2, 0.2, true); + ASSERT_EQ(elecstate->nelec_spin.size(), 2u); + EXPECT_NEAR(elecstate->nelec_spin[0], 4.2, 1e-14); + EXPECT_NEAR(elecstate->nelec_spin[1], 4.0, 1e-14); + EXPECT_TRUE(elecstate->eferm.two_efermi); + + elecstate::update_electron_count(*elecstate, 7.8, 2, -0.2, false); + EXPECT_NEAR(elecstate->nelec_spin[0], 3.8, 1e-14); + EXPECT_NEAR(elecstate->nelec_spin[1], 4.0, 1e-14); + EXPECT_FALSE(elecstate->eferm.two_efermi); + + elecstate::update_electron_count(*elecstate, 8.0, 2, 0.0, true); + EXPECT_DOUBLE_EQ(elecstate->nelec_spin[0], 4.0); + EXPECT_DOUBLE_EQ(elecstate->nelec_spin[1], 4.0); + EXPECT_TRUE(elecstate->eferm.two_efermi); + + // The explicit path must leave the native parameter source unchanged. + TestParameters::set_nspin(2); + elecstate->init_nelec_spin(); + EXPECT_DOUBLE_EQ(elecstate->nelec_spin[0], 5.0); + EXPECT_DOUBLE_EQ(elecstate->nelec_spin[1], 5.0); +} + TEST_F(ElecStateTest, Constructor) { Charge* charge = new Charge; diff --git a/source/source_estate/unittests/test_electronic_occupations.cpp b/source/source_estate/unittests/test_electronic_occupations.cpp new file mode 100644 index 00000000000..38fd39259bf --- /dev/null +++ b/source/source_estate/unittests/test_electronic_occupations.cpp @@ -0,0 +1,217 @@ +#include "../electronic_occupations.h" +#include "source_base/matrix.h" +#include "source_base/parallel_global.h" + +#include "gtest/gtest.h" + +#include +#include + +namespace +{ +int process_count = 1; +int process_rank = 0; + +ModuleBase::matrix fractional_weights() +{ + ModuleBase::matrix weights(2, 3); + weights(0, 0) = 2.0; + weights(0, 1) = 2.0; + weights(0, 2) = 2.0; + weights(1, 0) = 2.0; + weights(1, 1) = 0.2; + return weights; +} + +TEST(ElectronicOccupations, FractionalPopulationIncludesEveryLocalRowAndBand) +{ + const ModuleBase::matrix weights = fractional_weights(); + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, 1, process_count), 1e-14); +} + +TEST(ElectronicOccupations, ReturnsOccupationSumInsteadOfEchoingNearbyRequest) +{ + const ModuleBase::matrix weights = fractional_weights(); + const double requested = 8.2 + 5e-9; + const double actual = elecstate::checked_occupation_electron_count(weights, requested, 1, process_count); + EXPECT_NEAR(8.2, actual, 1e-14); + EXPECT_GT(requested - actual, 4e-9); +} + +TEST(ElectronicOccupations, RejectsPopulationMismatchOnEitherSide) +{ + const ModuleBase::matrix weights = fractional_weights(); + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2 + 2e-8, 1, process_count), + std::runtime_error); + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2 - 2e-8, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, ToleranceIsAbsoluteEvenForLargeElectronCounts) +{ + ModuleBase::matrix weights(1, 1); + weights(0, 0) = 1e6; + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 1e6 + 1e-6, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, GeneralizedSmearingCanContainNegativeWeights) +{ + ModuleBase::matrix weights(1, 3); + weights(0, 0) = 2.125; + weights(0, 1) = -0.25; + weights(0, 2) = 0.125; + EXPECT_DOUBLE_EQ(2.0, elecstate::checked_occupation_electron_count(weights, 2.0, 1, process_count)); +} + +TEST(ElectronicOccupations, RejectsNaNOnLastRankTogether) +{ + ModuleBase::matrix weights = fractional_weights(); + if (process_rank == process_count - 1) + { + weights(1, 2) = std::numeric_limits::quiet_NaN(); + } + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, RejectsInfinityOnLastRankTogether) +{ + ModuleBase::matrix weights = fractional_weights(); + if (process_rank == process_count - 1) + { + weights(1, 2) = std::numeric_limits::infinity(); + } + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, RejectsLocalOverflowBeforePoolReduction) +{ + ModuleBase::matrix weights(1, 2); + weights(0, 0) = std::numeric_limits::max(); + weights(0, 1) = std::numeric_limits::max(); + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, process_count, 1), + std::runtime_error); +} + +TEST(ElectronicOccupations, RejectsInvalidRequestOnLastRankTogether) +{ + const ModuleBase::matrix weights = fractional_weights(); + const double requested = process_rank == process_count - 1 ? std::numeric_limits::quiet_NaN() : 8.2; + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, requested, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, RejectsInvalidPartitionOnLastRankTogether) +{ + const ModuleBase::matrix weights = fractional_weights(); + const int local_size = process_rank == process_count - 1 ? 0 : process_count; + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, 1, local_size), std::runtime_error); + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, 0, process_count), std::runtime_error); +} + +TEST(ElectronicOccupations, SinglePoolDoesNotMultiplyReplicatedPopulationByRankCount) +{ + const ModuleBase::matrix weights = fractional_weights(); + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, 1, process_count), 1e-14); +} + +TEST(ElectronicOccupations, SinglePoolMismatchOnOneReplicaFailsEveryRank) +{ + ModuleBase::matrix weights = fractional_weights(); + if (process_rank == process_count - 1) + { + weights(1, 1) += 0.125; + } + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, 1, process_count), + std::runtime_error); +} + +TEST(ElectronicOccupations, KPointPoolsCombineDistinctRowsAndRemovePoolReplicas) +{ + const int pools = process_count % 2 == 0 ? 2 : 1; + const int ranks_per_pool = process_count / pools; + const int pool = process_rank / ranks_per_pool; + ModuleBase::matrix weights(1, 2); + const double population = pools == 1 ? 8.2 : (pool == 0 ? 2.2 : 6.0); + weights(0, 0) = 0.25 * population; + weights(0, 1) = 0.75 * population; + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, pools, ranks_per_pool), 1e-14); +} + +TEST(ElectronicOccupations, BandShardsUseLocalColumnCount) +{ + // Independent band shards are reduction groups, just like native BNDPAR. + const int local_bands = process_rank + 1; + ModuleBase::matrix weights(2, local_bands); + const double total_bands = process_count * (process_count + 1) / 2.0; + for (int band = 0; band < local_bands; ++band) + { + weights(0, band) = 3.0 / total_bands; + weights(1, band) = 5.2 / total_bands; + } + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, process_count, 1), 1e-14); +} + +TEST(ElectronicOccupations, EmptyLocalBlocksStillParticipateInCollectives) +{ + ModuleBase::matrix weights; + if (process_rank == 0) + { + weights = fractional_weights(); + } + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, process_count, 1), 1e-14); +} + +TEST(ElectronicOccupations, EmptyPopulationCannotSatisfyPositiveRequest) +{ + const ModuleBase::matrix weights; + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, process_count, 1), + std::runtime_error); +} + +TEST(ElectronicOccupations, UnequalPoolSizesRemoveTheirOwnReplicaCounts) +{ + // One rank in the first pool and all remaining ranks in the second. + const int pools = process_count > 1 ? 2 : 1; + const int ranks_per_pool = process_rank == 0 ? 1 : process_count - 1; + ModuleBase::matrix weights(1, 1); + weights(0, 0) = pools == 1 ? 8.2 : (process_rank == 0 ? 2.2 : 6.0); + EXPECT_NEAR(8.2, elecstate::checked_occupation_electron_count(weights, 8.2, pools, ranks_per_pool), 1e-14); +} + +TEST(ElectronicOccupations, RemoteInvalidEntryCannotBeHiddenByAnEmptyRootBlock) +{ + ModuleBase::matrix weights; + if (process_rank == process_count - 1) + { + weights.create(1, 1); + weights(0, 0) = -std::numeric_limits::infinity(); + } + EXPECT_THROW(elecstate::checked_occupation_electron_count(weights, 8.2, process_count, 1), + std::runtime_error); +} +} // namespace + +int main(int argc, char** argv) +{ + int thread_count = 1; + Parallel_Global::read_pal_param(argc, argv, process_count, thread_count, process_rank); +#ifdef __MPI + // Follow source_base/test/mpi_test_main.cpp: no pool communicators are + // created by these tests, so native finalization must not free them. + POOL_WORLD = MPI_COMM_NULL; + KP_WORLD = MPI_COMM_NULL; + INT_BGROUP = MPI_COMM_NULL; + BP_WORLD = MPI_COMM_NULL; + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + Parallel_Global::finalize_mpi(); +#endif + return result; +} diff --git a/source/source_estate/unittests/test_electrostatic_potential.cpp b/source/source_estate/unittests/test_electrostatic_potential.cpp new file mode 100644 index 00000000000..d7058c4f85f --- /dev/null +++ b/source/source_estate/unittests/test_electrostatic_potential.cpp @@ -0,0 +1,234 @@ +#include "source_estate/module_pot/electrostatic_potential.h" +#include "source_estate/module_pot/potential_new.h" +#include "source_estate/module_pot/efield.h" +#include "source_estate/module_pot/h_hartree_pw.h" +#include "source_base/parallel_global.h" +#include "source_io/module_parameter/input_parameter.h" + +#include +#include +#include +#include + +// Only the public cell/grid storage is needed by the averaging function. +UnitCell::UnitCell() {} +UnitCell::~UnitCell() {} +Magnetism::Magnetism() {} +Magnetism::~Magnetism() {} +SepPot::SepPot() {} +SepPot::~SepPot() {} +Sep_Cell::Sep_Cell() noexcept {} +Sep_Cell::~Sep_Cell() noexcept {} +Charge::Charge() {} +Charge::~Charge() {} +surchem::surchem() {} +surchem::~surchem() {} + +// A constant Hartree term isolates the assembly of the local potentials. +// The averaging tests supply their grid directly. +ModuleBase::matrix elecstate::H_Hartree_pw::v_hartree( + const UnitCell&, ModulePW::PW_Basis* basis, const int& nspin, const double* const* const) +{ + ModuleBase::matrix hartree(nspin, basis->nrxx); + for (int ir = 0; ir < basis->nrxx; ++ir) { hartree(0, ir) = 2.0; } + return hartree; +} +ModuleBase::matrix elecstate::Efield::add_efield( + const UnitCell&, const ModulePW::PW_Basis*, const int&, const double* const* const, const surchem&) +{ + throw std::logic_error("unexpected dipole construction in averaging test"); +} + +const std::vector* elecstate::Potential::pcc_electron_potential() const +{ + throw std::logic_error("unexpected PCC construction in averaging test"); +} + +const std::vector* elecstate::Potential::find_sol_pot() const +{ + throw std::logic_error("unexpected solvation construction in averaging test"); +} + +namespace +{ +int process_count = 1; +int process_rank = 0; + +class ElectrostaticPotentialTest : public testing::Test +{ + protected: + void set_grid(int planes) + { + cell.lat0_angstrom = 4.0; + const double z_length = planes / 4.0; + cell.latvec = ModuleBase::Matrix3(1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, z_length); + basis.nx = 4; + basis.ny = 4; + basis.nz = planes; + const int remainder = planes % process_count; + const int extra = process_rank < remainder ? 1 : 0; + basis.nplane = planes / process_count + extra; + basis.startz_current = process_rank * (planes / process_count) + std::min(process_rank, remainder); + basis.nrxx = basis.nx * basis.ny * basis.nplane; +#ifdef __MPI + basis.pool_world = MPI_COMM_WORLD; +#endif + potential.clear(); + for (int x = 0; x < basis.nx; ++x) + { + for (int y = 0; y < basis.ny; ++y) + { + for (int z = basis.startz_current; z < basis.startz_current + basis.nplane; ++z) + { + const double value = 100.0 * x + 10.0 * y + z; + potential.push_back(value); + } + } + } + } + + void expect_average(int axis, double lower, double upper, double expected) + { + const double actual = elecstate::planar_potential_average(cell, basis, potential, axis, lower, upper); + EXPECT_NEAR(actual, expected, 1.0e-12); + } + + UnitCell cell; + ModulePW::PW_Basis basis; + std::vector potential; +}; + +TEST_F(ElectrostaticPotentialTest, SccsAndPccAreAddedOnce) +{ + set_grid(5); + Input_para input; + input.nspin = 1; + input.imp_sol = 2; + Charge charge; + surchem solvent; + const std::vector fixed(basis.nrxx, 3.0); + const std::vector pcc(basis.nrxx, 0.25); + const std::vector sccs(basis.nrxx, 0.5); + // SCCS does not initialize the legacy delta_phi array. + ASSERT_EQ(solvent.delta_phi, nullptr); + const auto result = elecstate::electrostatic_potential(cell, basis, charge, fixed.data(), solvent, + input, &pcc, &sccs); + ASSERT_EQ(result.size(), static_cast(basis.nrxx)); + for (const double value : result) { EXPECT_DOUBLE_EQ(value, 5.75); } +} + +TEST_F(ElectrostaticPotentialTest, LegacySolventAndVacuumRemainUnchanged) +{ + set_grid(5); + Input_para input; + input.nspin = 1; + Charge charge; + surchem solvent; + const std::vector fixed(basis.nrxx, 3.0); + std::vector delta_phi(basis.nrxx, 0.75); + solvent.delta_phi = delta_phi.data(); + for (int model = 0; model <= 1; ++model) + { + input.imp_sol = model; + const auto result = elecstate::electrostatic_potential(cell, basis, charge, fixed.data(), solvent, + input, nullptr, nullptr); + const double expected = model == 0 ? 5.0 : 5.75; + for (const double value : result) { EXPECT_DOUBLE_EQ(value, expected); } + } +} + +TEST_F(ElectrostaticPotentialTest, SolventGridMismatchIsRejected) +{ + set_grid(5); + Input_para input; + input.nspin = 1; + input.imp_sol = 2; + Charge charge; + surchem solvent; + const std::vector fixed(basis.nrxx, 3.0); + const std::vector wrong_grid(basis.nrxx + 1, 0.5); + EXPECT_THROW(elecstate::electrostatic_potential(cell, basis, charge, fixed.data(), solvent, + input, nullptr, &wrong_grid), std::runtime_error); +} + +TEST_F(ElectrostaticPotentialTest, HalfOpenWindowsOnAllAxes) +{ + set_grid(5); + // V = 100*x + 10*y + z; unrestricted means are x=y=1.5, z=2. + expect_average(0, 1.0, 2.0, 117.0); + expect_average(1, 1.0, 2.0, 162.0); + expect_average(2, 1.0, 3.0, 166.5); +} + +TEST_F(ElectrostaticPotentialTest, UnequalPartitionsUsePointCounts) +{ + set_grid(5); + // Four ranks hold 2/1/1/1 planes; averaging rank means would be wrong. + expect_average(2, 0.0, 5.0, 167.0); +} + +TEST_F(ElectrostaticPotentialTest, EmptyGlobalWindowFailsOnEveryRank) +{ + set_grid(5); + EXPECT_THROW(elecstate::planar_potential_average(cell, basis, potential, 2, 0.1, 0.2), + std::runtime_error); +} + +TEST_F(ElectrostaticPotentialTest, EmptyLocalPlanesContributeZero) +{ + set_grid(2); + // With four ranks, two ranks own no planes. + if (process_rank >= 2) + { + EXPECT_EQ(basis.nplane, 0); + EXPECT_TRUE(potential.empty()); + } + expect_average(2, 0.0, 2.0, 165.5); +} + +TEST_F(ElectrostaticPotentialTest, RanksOutsideWindowStillParticipate) +{ + set_grid(5); + expect_average(2, 4.0, 5.0, 169.0); +} + +TEST_F(ElectrostaticPotentialTest, TriclinicWindowsUseNormalHeights) +{ + set_grid(5); + cell.latvec = ModuleBase::Matrix3(1.0, 0.0, 0.0, + 0.75, 1.0, 0.0, + 0.75, 0.75, 1.0); + // For this cell the normal heights are 64/sqrt(409), 3.2, and 4 A. + const double heights[3] = {64.0 / std::sqrt(409.0), 3.2, 4.0}; + const double centers[3] = {0.25, 0.25, 0.2}; + const double expected[3] = {117.0, 162.0, 166.0}; + for (int axis = 0; axis < 3; ++axis) + { + const double lower = (centers[axis] - 0.01) * heights[axis]; + const double upper = (centers[axis] + 0.01) * heights[axis]; + expect_average(axis, lower, upper, expected[axis]); + } +} +} // namespace + +int main(int argc, char** argv) +{ + int thread_count = 1; + Parallel_Global::read_pal_param(argc, argv, process_count, thread_count, process_rank); +#ifdef __MPI + POOL_WORLD = MPI_COMM_NULL; + KP_WORLD = MPI_COMM_NULL; + INT_BGROUP = MPI_COMM_NULL; + BP_WORLD = MPI_COMM_NULL; + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + Parallel_Global::finalize_mpi(); +#endif + return result; +} diff --git a/source/source_estate/write_elecstat_pot.cpp b/source/source_estate/write_elecstat_pot.cpp index f530c3d3d0c..68fc265cda4 100644 --- a/source/source_estate/write_elecstat_pot.cpp +++ b/source/source_estate/write_elecstat_pot.cpp @@ -1,8 +1,7 @@ #include "source_base/element_name.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" -#include "source_estate/module_pot/h_hartree_pw.h" -#include "source_estate/module_pot/efield.h" +#include "source_estate/module_pot/electrostatic_potential.h" #include "source_io/module_output/cube_io.h" #include "source_io/module_output/output_log.h" #include "write_elecstat_pot.h" @@ -30,63 +29,9 @@ void write_elecstat_pot( ModuleBase::TITLE("ModuleIO", "write_elecstat_pot"); ModuleBase::timer::start("ModuleIO", "write_elecstat_pot"); - std::vector v_elecstat(rho_basis->nrxx, 0.0); - - if (solvent_potential != nullptr && solvent_potential->size() != static_cast(rho_basis->nrxx)) - { - ModuleBase::WARNING_QUIT("write_elecstat_pot", "Solvent potential does not match the output grid"); - } - - const int nspin = PARAM.inp.nspin; - const int efield = PARAM.inp.efield_flag; - const int dip_corr = PARAM.inp.dip_cor_flag; - const int imp_sol = PARAM.inp.imp_sol; - - //========================================== - // Hartree potential - //========================================== - ModuleBase::matrix vh(nspin, rho_basis->nrxx); - vh = elecstate::H_Hartree_pw::v_hartree(*ucell, rho_basis, nspin, chr->rho); - - //========================================== - //! Dipole correction - //========================================== - ModuleBase::matrix v_efield; - if (efield>0 && dip_corr>0) - { - v_efield.create(nspin, rho_basis->nrxx); - v_efield = elecstate::Efield::add_efield(*ucell, - const_cast(rho_basis), - nspin, - chr->rho, - solvent); - } - - //========================================== - //! Add hartree potential and local pseudopot - //========================================== - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - // the spin index is 0 - v_elecstat[ir] = vh(0, ir) + v_eff[ir]; - - if (efield>0 && dip_corr>0) - { - v_elecstat[ir] += v_efield(0, ir); - } - if (pcc_potential != nullptr) - { - v_elecstat[ir] += (*pcc_potential)[ir]; - } - if (solvent_potential != nullptr) - { - v_elecstat[ir] += (*solvent_potential)[ir]; - } - if (imp_sol == 1) - { - v_elecstat[ir] += solvent.delta_phi[ir]; - } - } + auto v_elecstat = elecstate::electrostatic_potential(*ucell, *rho_basis, *chr, + v_eff, solvent, PARAM.inp, + pcc_potential, solvent_potential); //------------------------------------------- //! Get the vacuum level of the system @@ -113,7 +58,7 @@ void write_elecstat_pot( ModuleIO::write_vdata_palgrid(*chr->pgrid, v_elecstat.data(), is, - nspin, + PARAM.inp.nspin, istep, fn, ef_tmp, diff --git a/source/source_hsolver/unittests/hsolver_supplementary_mock.h b/source/source_hsolver/unittests/hsolver_supplementary_mock.h index 6d216e7ec81..3b10abc1e72 100644 --- a/source/source_hsolver/unittests/hsolver_supplementary_mock.h +++ b/source/source_hsolver/unittests/hsolver_supplementary_mock.h @@ -15,6 +15,8 @@ const double *ElecState::getRho(int spin) const void ElecState::init_nelec_spin() { return; } +void ElecState::init_nelec_spin(int, double, double) {} + void ElecState::init_ks( Charge *chg_in, // pointer for class Charge const K_Vectors *klist_in, diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index 51c8d95a8b6..791ce76e7e9 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -222,7 +222,11 @@ void ReadInput::item_system() 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.)"; +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 ordinary SCF step is returned with scf_converged=false metadata so an external driver can choose its policy. + +Socket mode defaults symmetry to 0 so moving atoms do not retain the initial STRU spatial symmetry; an explicit symmetry input remains unchanged. + +Constant-potential i-PI sessions use INIT electron metadata and CHGDATA before each POSDATA to set the total electron number. Use fixed-cell CPU double-precision norm-conserving PW/LCAO KSDFT with cal_force=1, cal_stress=0, symmetry=0 or -1, and nspin=1 or 2. CHGDATA flags are zero. Constant-potential electron control is incompatible with the energy-only Makov-Payne correction (assume_isolated=makov-payne). Each response contains the converged energy and forces, actual electron count, Fermi level, and requested work function. Set electron-control options in the i-PI electrons block. See the [i-PI examples and tests](https://github.com/zhangchenyu-USTC/i-pi_constant-potential/tree/abacus-develop-validation/examples/clients/abacus/validation).)"; item.default_value = "False"; read_sync_bool(input.socket_driver); item.check_value = [](const Input_Item& item, const Parameter& para) { @@ -340,7 +344,7 @@ Socket mode always computes energy. Force and stress extraction follows cal_forc { para.input.symmetry = "-1"; } - else if (para.input.gamma_only || para.input.calculation == "nscf" || para.input.calculation == "get_s" + else if (para.input.socket_driver || para.input.gamma_only || para.input.calculation == "nscf" || para.input.calculation == "get_s" || para.input.calculation == "get_pchg" || para.input.calculation == "get_wf") { para.input.symmetry = "0"; // if md or exx, symmetry will be 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 e66274bc7c8..1aeb96f4e80 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -375,6 +375,14 @@ TEST_F(InputTest, Item_test) } { // symmetry auto it = find_label("symmetry", readinput.input_lists); + TestParameters::input(param).socket_driver = true; + TestParameters::input(param).symmetry = "default"; + it->second.reset_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).symmetry, "0"); + TestParameters::input(param).symmetry = "1"; + it->second.reset_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).symmetry, "1"); + TestParameters::input(param).socket_driver = false; TestParameters::input(param).symmetry = "default"; TestParameters::input(param).gamma_only = true; it->second.reset_value(it->second, param); diff --git a/source/source_lcao/module_ri/exx_lri_interface.h b/source/source_lcao/module_ri/exx_lri_interface.h index 8b5bc61cf64..9a55f16b1e7 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.h +++ b/source/source_lcao/module_ri/exx_lri_interface.h @@ -99,6 +99,9 @@ class Exx_LRI_Interface void exx_beforescf(const int istep, const K_Vectors& kv, const Charge_Mixing& chgmix, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_2D& pv, const bool update_symmetry); + // Only the explicit electronic-control entry point enables strict EXX convergence. + void begin_electronic_step() { strict_electronic_step_ = true; } + /// @brief in eachiterinit: do DM mixing and calculate Hexx when entering 2nd SCF void exx_eachiterinit(const int istep, const UnitCell& ucell, @@ -139,6 +142,7 @@ class Exx_LRI_Interface /// 0: the first outer loop. If restart, contributeHk will directly add Hexx to Hloc. else, do nothing. /// >0: not the first outer loop. contributeHk will do enerything normally. int two_level_step = 0; + bool strict_electronic_step_ = false; double etot_last_outer_loop = 0.0; module_dm::DensityMatrix* dm_last_step; diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index b08514faad0..ca279b1981d 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -411,6 +411,10 @@ bool Exx_LRI_Interface::exx_after_converge( const double ediff = std::abs(etot - etot_last_outer_loop) * ModuleBase::Ry_to_eV; if (two_level_step) { std::cout << FmtCore::format(" deltaE (eV) from outer loop: %.8e \n", ediff); } + const bool outer_converged = this->two_level_step != 0 + && (iter == 1 || ediff < scf_ene_thr); + if (strict_electronic_step_ && this->two_level_step >= this->hybrid_step_ && !outer_converged) + { throw std::runtime_error("constant-potential EXX outer loop did not converge; increase exx_hybrid_step"); } // exx converged or get max exx steps if (this->two_level_step == this->hybrid_step_ || (iter == 1 && this->two_level_step != 0) // density convergence of outer loop diff --git a/source/source_main/driver_run.cpp b/source/source_main/driver_run.cpp index ccfa6e55959..69fa638aabc 100644 --- a/source/source_main/driver_run.cpp +++ b/source/source_main/driver_run.cpp @@ -114,7 +114,7 @@ void Driver::driver_run() 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); + rl_driver.relax_driver(p_esolver, ucell, PARAM, GlobalV::ofs_running); } else if (cal == "get_s") { diff --git a/source/source_pw/module_pwdft/exx_helper.cpp b/source/source_pw/module_pwdft/exx_helper.cpp index 79d102becd1..45fd262ad74 100644 --- a/source/source_pw/module_pwdft/exx_helper.cpp +++ b/source/source_pw/module_pwdft/exx_helper.cpp @@ -7,6 +7,19 @@ #include "source_estate/elecstate_pw.h" // use ElecStatePW #include "source_estate/module_charge/charge.h" // use Charge #include // for timing +#include + +void validate_electronic_step_exx(int functional_type, bool cal_exx) +{ + // Separate-loop SCAN0 starts as SCAN, before switching to hybrid meta-GGA. + // The upstream PW density path updates tau only for non-hybrid meta-GGA. + if (functional_type == 5 || (functional_type == 3 && cal_exx)) + { + throw std::runtime_error( + "constant-potential PW hybrid meta-GGA is unsupported: " + "upstream kinetic-energy-density updates do not handle hybrid functionals"); + } +} template void Exx_Helper::init(const UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& wg, const General_Exx_Info& exx_info) @@ -90,7 +103,7 @@ bool Exx_Helper::iter_finish(void* p_elec, Charge* p_charge, void* ps } bool conv_ene = std::abs(dexx) < inp.exx_ene_thr; - conv_esolver = this->exx_after_converge(iter, conv_ene); + conv_esolver = this->exx_after_converge(iter, conv_ene, inp); if (!conv_esolver) { @@ -130,9 +143,20 @@ double Exx_Helper::cal_exx_energy(void* psi_) } template -bool Exx_Helper::exx_after_converge(int iter, bool ene_conv) +bool Exx_Helper::exx_after_converge(int iter, bool ene_conv, const Input_para& inp) { - if (op_exx->first_iter) + if (strict_electronic_step && op_exx->separate_loop) + { + // A changed electron number must complete at least one EXX update. + // Reaching the outer-loop budget is not a convergence criterion. + const bool outer_converged = exx_iter > 0 + && ((inp.exx_thr_type == "energy" && ene_conv) + || (inp.exx_thr_type == "density" && iter == 1)); + if (outer_converged) { return true; } + if (exx_iter >= inp.exx_hybrid_step) + { throw std::runtime_error("constant-potential EXX outer loop did not converge; increase exx_hybrid_step"); } + } + else if (op_exx->first_iter) { op_exx->first_iter = false; } @@ -140,15 +164,15 @@ bool Exx_Helper::exx_after_converge(int iter, bool ene_conv) { return true; } - else if (PARAM.inp.exx_thr_type == "energy" && ene_conv) + else if (inp.exx_thr_type == "energy" && ene_conv) { return true; } - else if (PARAM.inp.exx_thr_type == "density" && iter == 1) + else if (inp.exx_thr_type == "density" && iter == 1) { return true; } - else if (exx_iter >= PARAM.inp.exx_hybrid_step) + else if (exx_iter >= inp.exx_hybrid_step) { GlobalV::ofs_running << " !!EXX IS NOT CONVERGED!!" << std::endl; std::cout << " !!EXX IS NOT CONVERGED!!" << std::endl; @@ -188,4 +212,3 @@ template class Exx_Helper, base_device::DEVICE_CPU>; template class Exx_Helper, base_device::DEVICE_GPU>; template class Exx_Helper, base_device::DEVICE_GPU>; #endif - diff --git a/source/source_pw/module_pwdft/exx_helper.h b/source/source_pw/module_pwdft/exx_helper.h index 31b25eeeadf..8f7fb3d8ebb 100644 --- a/source/source_pw/module_pwdft/exx_helper.h +++ b/source/source_pw/module_pwdft/exx_helper.h @@ -20,6 +20,12 @@ struct Exx_Helper : public Exx_HelperBase virtual ~Exx_Helper() = default; OperatorEXX *op_exx = nullptr; + void begin_electronic_step() override + { + strict_electronic_step = true; + exx_iter = 0; + } + void init(const UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& wg, const General_Exx_Info& exx_info) override; void before_scf(void* p_hamilt, void* psi, const Input_para& inp, const General_Exx_Info& exx_info) override; @@ -36,11 +42,16 @@ struct Exx_Helper : public Exx_HelperBase void set_op() override { op_exx->first_iter = first_iter; + if (strict_electronic_step) + { + // Warm starts construct ACE immediately and require valid occupations. + op_exx->set_wg(wg); + } set_psi(psi); - op_exx->set_wg(wg); + if (!strict_electronic_step) { op_exx->set_wg(wg); } } - bool exx_after_converge(int iter, bool ene_conv) override; + bool exx_after_converge(int iter, bool ene_conv, const Input_para& inp) override; double cal_exx_energy(void* psi_) override; @@ -53,6 +64,7 @@ struct Exx_Helper : public Exx_HelperBase psi::Psi *psi = nullptr; const ModuleBase::matrix *wg = nullptr; int exx_iter = 0; + bool strict_electronic_step = false; }; #endif // EXX_HELPER_H diff --git a/source/source_pw/module_pwdft/exx_helper_base.h b/source/source_pw/module_pwdft/exx_helper_base.h index 12bd8d7aeea..11f1b98b53a 100644 --- a/source/source_pw/module_pwdft/exx_helper_base.h +++ b/source/source_pw/module_pwdft/exx_helper_base.h @@ -8,12 +8,17 @@ class UnitCell; struct Input_para; struct General_Exx_Info; +void validate_electronic_step_exx(int functional_type, bool cal_exx); + class Exx_HelperBase { public: Exx_HelperBase() = default; virtual ~Exx_HelperBase() = default; + // Explicit boundary for an externally controlled electron-number request. + virtual void begin_electronic_step() = 0; + virtual void init(const UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& wg, const General_Exx_Info& exx_info) = 0; virtual void before_scf(void* p_hamilt, void* psi, const Input_para& inp, const General_Exx_Info& exx_info) = 0; @@ -35,7 +40,7 @@ class Exx_HelperBase /// value (false = not converged, rerun SCF) and the caller /// (ESolver_KS_PW) sets ESolver_KS::scf_rerun_ to trigger the /// restart at iter=1 in ESolver_KS::runner. - virtual bool exx_after_converge(int iter, bool ene_conv) = 0; + virtual bool exx_after_converge(int iter, bool ene_conv, const Input_para& inp) = 0; virtual double cal_exx_energy(void* psi) = 0; diff --git a/source/source_relax/CMakeLists.txt b/source/source_relax/CMakeLists.txt index 8dc901592f7..0f0f38f7514 100644 --- a/source/source_relax/CMakeLists.txt +++ b/source/source_relax/CMakeLists.txt @@ -6,6 +6,7 @@ add_library( socket_frame.cpp socket_utils.cpp socket_handlers.cpp + socket_constant_potential.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 cf56ffd2da4..096a9dd2102 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -12,10 +12,12 @@ #include "source_io/module_parameter/parameter.h" void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - std::ofstream& ofs_running) + UnitCell& ucell, + Parameter& parameters, + std::ofstream& ofs_running) { + const Input_para& inp = parameters.inp; + const System_para& system = parameters.globalv; ModuleBase::TITLE("Relax_Driver", "relax_driver"); ModuleBase::timer::start("Relax_Driver", "relax_driver"); @@ -28,8 +30,8 @@ void Relax_Driver::relax_driver(ModuleESolver::ESolver* p_esolver, if (inp.socket_driver) { Socket_Driver socket_driver; - socket_driver.socket_driver(p_esolver, ucell, inp, ofs_running); - p_esolver->after_all_runners(ucell); + const bool ran_frame = socket_driver.socket_driver(p_esolver, ucell, parameters, inp, system, ofs_running); + if (ran_frame) { p_esolver->after_all_runners(ucell); } ModuleBase::timer::end("Relax_Driver", "relax_driver"); return; } diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index 34540f00b7f..2f02ac11412 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -18,6 +18,8 @@ * This class manages the main relaxation loop, including energy solving, * force/stress calculation, and structure optimization steps. */ +class Parameter; + class Relax_Driver { @@ -34,12 +36,12 @@ class Relax_Driver * * @param p_esolver Pointer to the energy solver. * @param ucell Reference to the unit cell to be relaxed. - * @param inp Input parameters for the calculation. + * @param parameters Runtime input and electronic parameters. * @param ofs_running Output stream for running log. */ void relax_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell, - const Input_para& inp, + Parameter& parameters, std::ofstream& ofs_running); private: diff --git a/source/source_relax/socket_constant_potential.cpp b/source/source_relax/socket_constant_potential.cpp new file mode 100644 index 00000000000..cecff4c7d0d --- /dev/null +++ b/source/source_relax/socket_constant_potential.cpp @@ -0,0 +1,198 @@ +#include "socket_constant_potential.h" +#include "source_base/constants.h" +#include "source_esolver/esolver.h" +#include "source_io/module_parameter/input_parameter.h" +#include "source_io/module_parameter/system_parameter.h" + +#include +#include +#include +#include +#include +#include + +namespace +{ +std::string trim(const std::string& text) +{ + const auto begin = text.find_first_not_of(" \t\r\n"); + if (begin == std::string::npos) { return ""; } + return text.substr(begin, text.find_last_not_of(" \t\r\n") - begin + 1); +} + +double number(const std::map& tokens, const std::string& key) +{ + const auto found = tokens.find(key); + if (found == tokens.end()) { throw std::runtime_error("constant-potential INIT requires " + key); } + std::istringstream in(found->second); + in.imbue(std::locale::classic()); + double value = 0.0; + if (!(in >> value) || !std::isfinite(value)) + { throw std::runtime_error("INIT " + key + " must be finite and numerical"); } + in >> std::ws; + if (!in.eof()) { throw std::runtime_error("invalid trailing data in INIT " + key); } + return value; +} +} + +namespace SocketCP +{ +void run_electronic_step(ModuleESolver::ESolver& solver, UnitCell& cell, + Parameter& parameters, const Input_para& inp, const System_para& system, + const elecstate::ElectronicControl& control, int functional_type, int istep) +{ + if (inp.assume_isolated == "makov-payne") + { + throw std::runtime_error( + "constant-potential electron control is incompatible with the energy-only Makov-Payne correction"); + } + if (inp.esolver_type != "ksdft" || (inp.basis_type != "pw" && inp.basis_type != "lcao") + || inp.device != "cpu" || inp.precision != "double" || system.use_uspp + || inp.sc_mag_switch || inp.deepks_scf || inp.rpa || inp.rdmft || inp.dft_plus_dmft) + { + throw std::runtime_error( + "constant-potential control requires CPU double-precision PW/LCAO KSDFT with " + "norm-conserving pseudopotentials; specialized solver corrections are not supported"); + } + if (functional_type < 1 || functional_type > 5) + { + throw std::runtime_error( + "constant-potential control requires a supported KS exchange-correlation functional"); + } + if (inp.dft_plus_u != 0 && inp.dft_plus_u != 1) + { + throw std::runtime_error( + "constant-potential forces require dft_plus_u=1; " + "the legacy dft_plus_u=2 force implementation is unavailable"); + } + if (inp.dft_plus_u && (inp.init_occ_mat < 0 || inp.init_occ_mat > 1 + || inp.yukawa_potential || system.uramping > 0.0)) + { + throw std::runtime_error( + "constant-potential DFT+U requires fixed Hubbard U and self-consistent occupations " + "(init_occ_mat=0/1, yukawa_potential=false, uramping disabled)"); + } + if (inp.init_chg == "dm_no_renormalize") + { throw std::runtime_error("electron-number control requires charge-density renormalization"); } + if (inp.nspin != 1 && inp.nspin != 2) + { throw std::runtime_error("constant-potential control supports only nspin=1 or 2"); } + if (inp.ocp || inp.nbands <= 0 || !std::isfinite(control.nelect) || control.nelect <= 0.0) + { throw std::runtime_error("invalid total electron number, band capacity, or fixed occupation array"); } + const bool two_fermi = system.two_fermi; + if (two_fermi && inp.nspin != 2) + { throw std::runtime_error("two_fermi electron control requires nspin=2"); } + const double up = 0.5 * (control.nelect + inp.nupdown); + const double down = 0.5 * (control.nelect - inp.nupdown); + if ((two_fermi && (!std::isfinite(inp.nupdown) || up <= 0.0 || down <= 0.0 + || up > inp.nbands || down > inp.nbands)) + || control.nelect > 2.0 * inp.nbands) + { + throw std::runtime_error( + "electron number exceeds band capacity or gives invalid spin populations; increase nbands"); + } + if (inp.smearing_method == "fixed" + && (std::abs(control.nelect - std::round(control.nelect)) > 1.0e-8 + || (two_fermi && (std::abs(up-std::round(up)) > 1.0e-8 || std::abs(down-std::round(down)) > 1.0e-8)))) + { throw std::runtime_error("fractional electron populations require a smearing occupation method"); } + if (!solver.apply_electronic_control(parameters, control)) + { + throw std::runtime_error("selected ESolver does not support constant-potential electron control"); + } + solver.runner(cell, istep); +} + +Configuration parse_init(const std::string& parameters) +{ + const std::string keys[] = {"NELECT", "neutral_electrons", + "potential_average_axis", "potential_average_min_A", "potential_average_max_A"}; + std::map tokens; + std::istringstream input(parameters); + std::string token; + while (std::getline(input, token, ',')) + { + token = trim(token); + const auto colon = token.find(':'); + const std::string key = trim(token.substr(0, colon)); + if (key == "solvation_update_stride") + { throw std::runtime_error("INIT solvation_update_stride has been removed"); } + for (const auto& known : keys) + { + if (key != known) { continue; } + if (colon == std::string::npos || tokens.count(key)) + { throw std::runtime_error("missing separator or duplicate constant-potential INIT token: " + key); } + tokens[key] = trim(token.substr(colon + 1)); + } + } + Configuration config; + if (tokens.empty()) { return config; } + config.enabled = true; + if (number(tokens, "NELECT") <= 0.0) + { throw std::runtime_error("INIT NELECT must be positive"); } + if (number(tokens, "neutral_electrons") <= 0.0) + { throw std::runtime_error("INIT neutral_electrons must be positive"); } + if (tokens.count("potential_average_axis") || tokens.count("potential_average_min_A") + || tokens.count("potential_average_max_A")) + { + const std::string axis = tokens["potential_average_axis"]; + if (axis.size() != 1 || axis[0] < 'x' || axis[0] > 'z') + { throw std::runtime_error("INIT potential_average_axis must be x, y or z"); } + config.window.axis = axis[0] - 'x'; + config.window.minimum_A = number(tokens, "potential_average_min_A"); + config.window.maximum_A = number(tokens, "potential_average_max_A"); + if (config.window.minimum_A < 0.0 || config.window.maximum_A <= config.window.minimum_A) + { throw std::runtime_error("INIT potential window must satisfy 0 <= min < max"); } + } + return config; +} + +elecstate::ElectronicControl make_control(const Configuration& config, double nelect, int flags) +{ + if (!config.enabled) { throw std::runtime_error("CHGDATA requires constant-potential INIT metadata"); } + if (!std::isfinite(nelect) || nelect <= 0.0) + { throw std::runtime_error("CHGDATA total electron number must be finite and positive"); } + if (flags != 0 && flags != 2) + { throw std::runtime_error("CHGDATA flags must be 0 or 2; bit 0 is reserved"); } + if (flags & 2) { throw std::runtime_error("ABACUS Ne-electrode CHGDATA is not supported"); } + elecstate::ElectronicControl control; + control.nelect = nelect; + return control; +} + +Observables collect_observables(const Configuration& config, const ModuleESolver::ESolver& solver, + const UnitCell& cell, const Parameter& parameters) +{ + if (!config.enabled) { throw std::runtime_error("electronic observables require constant-potential INIT"); } + elecstate::ElectronicObservables physical; + if (!solver.electronic_observables(cell, parameters, config.window, physical)) + { throw std::runtime_error("ESolver cannot provide constant-potential electronic observables"); } + Observables result; + result.fermi_eV = physical.fermi_ry * ModuleBase::Ry_to_eV; + result.nelect = physical.nelect; + if (!std::isfinite(result.fermi_eV) || !std::isfinite(result.nelect) || result.nelect <= 0.0) + { throw std::runtime_error("invalid constant-potential electronic observables"); } + result.has_workfunction = config.window.axis >= 0; + if (result.has_workfunction) + { + result.reference_eV = physical.reference_ry * ModuleBase::Ry_to_eV; + result.workfunction_eV = result.reference_eV - result.fermi_eV; + if (!std::isfinite(result.reference_eV) || !std::isfinite(result.workfunction_eV)) + { throw std::runtime_error("non-finite workfunction"); } + } + return result; +} + +std::string serialize_observables(const Observables& result) +{ + std::ostringstream out; + out.imbue(std::locale::classic()); + out << std::setprecision(17) << "{\"fermi_level_eV\":" << result.fermi_eV + << ",\"nelect\":" << result.nelect; + if (result.has_workfunction) + { + out << ",\"workfunction_eV\":" << result.workfunction_eV + << ",\"electrostatic_reference_eV\":" << result.reference_eV; + } + out << "}"; + return out.str(); +} +} diff --git a/source/source_relax/socket_constant_potential.h b/source/source_relax/socket_constant_potential.h new file mode 100644 index 00000000000..7ffc1d9483e --- /dev/null +++ b/source/source_relax/socket_constant_potential.h @@ -0,0 +1,44 @@ +#ifndef ABACUS_SOCKET_CONSTANT_POTENTIAL_H +#define ABACUS_SOCKET_CONSTANT_POTENTIAL_H + +#include "source_estate/electronic_control.h" +#include + +class Parameter; +class UnitCell; +struct Input_para; +struct System_para; +namespace ModuleESolver { class ESolver; } + +namespace SocketCP +{ +struct Configuration +{ + bool enabled = false; + elecstate::PotentialWindow window; +}; + +struct Observables +{ + double fermi_eV = 0.0; + double nelect = 0.0; + bool has_workfunction = false; + double reference_eV = 0.0; + double workfunction_eV = 0.0; +}; + +// Apply one external electron-number request and evaluate the current geometry. +void run_electronic_step(ModuleESolver::ESolver& solver, UnitCell& cell, + Parameter& parameters, const Input_para& input, const System_para& system, + const elecstate::ElectronicControl& control, int functional_type, int istep); + +Configuration parse_init(const std::string& parameters); +elecstate::ElectronicControl make_control(const Configuration& config, double nelect, int flags); +// Called on every rank while evaluating POSDATA; the solver may use collectives. +Observables collect_observables(const Configuration& config, const ModuleESolver::ESolver& solver, + const UnitCell& cell, const Parameter& parameters); +// Serializes only the completed frame, without calling the solver or MPI. +std::string serialize_observables(const Observables& result); +} + +#endif diff --git a/source/source_relax/socket_driver.cpp b/source/source_relax/socket_driver.cpp index 03d19cb2764..433c2688e3d 100644 --- a/source/source_relax/socket_driver.cpp +++ b/source/source_relax/socket_driver.cpp @@ -51,10 +51,12 @@ void log_peer_closed(std::ofstream& ofs_running) } } // namespace -void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, - UnitCell& ucell, - const Input_para& inp, - std::ofstream& ofs_running) +static bool run_socket_driver(ModuleESolver::ESolver* p_esolver, + UnitCell& ucell, + const Input_para& inp, + Parameter* parameters, + const System_para* system, + std::ofstream& ofs_running) { ModuleBase::TITLE("Socket_Driver", "socket_driver"); ModuleBase::timer::start("Socket_Driver", "socket_driver"); @@ -64,22 +66,24 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, ModuleBase::WARNING_QUIT("ABACUS socket", "socket driver requires a valid ESolver."); } IpiSocket socket; + DriverContext context; try { connect_on_root(socket, ofs_running); - DriverContext context; context.esolver = p_esolver; context.ucell = &ucell; context.inp = &inp; + context.parameters = parameters; + context.system = system; context.state = DriverState::NeedInit; context.nat_return = ucell.nat; context.reference_cell = ipi_cell_bohr_from_unitcell(ucell); while (true) { - const std::string header = SocketHandlers::read_header_bcast(socket, context.state); + const std::string header = SocketHandlers::read_header_bcast(socket, context); if (header.empty()) { log_peer_closed(ofs_running); @@ -97,13 +101,17 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, { SocketHandlers::handle_posdata(socket, context, ofs_running); } + else if (header == "CHGDATA") + { + SocketHandlers::handle_chgdata(socket, context, ofs_running); + } else if (header == "GETFORCE") { - SocketHandlers::handle_getforce(socket, context); + SocketHandlers::handle_getforce(socket, context, ofs_running); } else if (header == "EXIT") { - SocketHandlers::handle_exit(ofs_running); + SocketHandlers::handle_exit(context, ofs_running); break; } else @@ -123,4 +131,18 @@ void Socket_Driver::socket_driver(ModuleESolver::ESolver* p_esolver, } ModuleBase::timer::end("Socket_Driver", "socket_driver"); + return context.istep > 0; +} + +bool Socket_Driver::socket_driver(ModuleESolver::ESolver* solver, UnitCell& cell, + const Input_para& input, std::ofstream& log) +{ + return run_socket_driver(solver, cell, input, nullptr, nullptr, log); +} + +bool Socket_Driver::socket_driver(ModuleESolver::ESolver* solver, UnitCell& cell, + Parameter& parameters, const Input_para& input, + const System_para& system, std::ofstream& log) +{ + return run_socket_driver(solver, cell, input, ¶meters, &system, log); } diff --git a/source/source_relax/socket_driver.h b/source/source_relax/socket_driver.h index 86c180fc42e..8a37dc44039 100644 --- a/source/source_relax/socket_driver.h +++ b/source/source_relax/socket_driver.h @@ -4,7 +4,9 @@ #include class UnitCell; +class Parameter; struct Input_para; +struct System_para; namespace ModuleESolver { @@ -17,7 +19,11 @@ class Socket_Driver Socket_Driver() = default; ~Socket_Driver() = default; - void socket_driver(ModuleESolver::ESolver* p_esolver, + bool socket_driver(ModuleESolver::ESolver* solver, UnitCell& cell, + Parameter& parameters, const Input_para& input, + const System_para& system, std::ofstream& log); + + bool socket_driver(ModuleESolver::ESolver* p_esolver, UnitCell& ucell, const Input_para& inp, std::ofstream& ofs_running); diff --git a/source/source_relax/socket_handlers.cpp b/source/source_relax/socket_handlers.cpp index 71ceb4c4eab..e89e87fc133 100644 --- a/source/source_relax/socket_handlers.cpp +++ b/source/source_relax/socket_handlers.cpp @@ -3,6 +3,7 @@ #include "source_base/global_function.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" +#include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/input_parameter.h" #include "source_relax/socket_frame.h" #ifdef __JSON @@ -12,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -36,6 +38,7 @@ using SocketUtils::kMaxInitBytes; using SocketUtils::kRyToHartree; using SocketUtils::kStressAbsoluteTolerance; using SocketUtils::kStressRelativeTolerance; +using SocketUtils::bcast_string; using SocketUtils::matrix9_from_stress; using SocketUtils::max_abs_delta; using SocketUtils::max_wrapped_direct_delta_from_unitcell; @@ -46,7 +49,7 @@ 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 read_header_bcast(IpiSocket& socket, const DriverContext& context) { std::string header; int io_failed = 0; @@ -59,7 +62,7 @@ std::string read_header_bcast(IpiSocket& socket, const DriverState state) } catch (const IpiSocketClosed&) { - if (state == DriverState::HasData) + if (context.state == DriverState::HasData || context.charge_pending) { io_failed = 1; io_message = "i-PI peer closed while a computed frame was pending"; @@ -154,6 +157,25 @@ void handle_init(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_r { ofs_running << " ABACUS socket INIT params bytes " << nbytes << std::endl; } + bcast_string(params); + int local_failed = 0; + std::string local_message; + try + { + context.cp = SocketCP::parse_init(params); + if (context.cp.enabled && (context.parameters == nullptr || context.system == nullptr || !context.inp->cal_force + || context.inp->cal_stress || (context.inp->symmetry != "0" && context.inp->symmetry != "-1"))) + { + throw std::runtime_error("constant-potential socket requires runtime parameters, " + "cal_force=1, cal_stress=0 and symmetry=0/-1"); + } + } + catch (const std::exception& exc) + { + local_failed = 1; + local_message = exc.what(); + } + throw_if_any_rank_failed(local_failed, local_message); context.state = DriverState::Ready; if (is_root()) { @@ -161,6 +183,52 @@ void handle_init(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_r } } +void handle_chgdata(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running) +{ + std::vector charge(1, 0.0); + std::int32_t flags = 0; + int io_failed = 0; + std::string io_message; + if (is_root()) + { + try + { + if (context.state != DriverState::Ready || context.charge_pending) + { + throw std::runtime_error("CHGDATA requires READY with no pending charge or force result"); + } + charge[0] = socket.read_double(); + flags = socket.read_int32(); + } + catch (const std::exception& exc) + { + io_failed = 1; + io_message = exc.what(); + } + } + quit_if_root_failed(io_failed, io_message); + bcast_double_vector(charge); + bcast_int32(flags); + int local_failed = 0; + std::string local_message; + try + { + context.pending_control = SocketCP::make_control(context.cp, charge[0], flags); + } + catch (const std::exception& exc) + { + local_failed = 1; + local_message = exc.what(); + } + throw_if_any_rank_failed(local_failed, local_message); + context.charge_pending = true; + if (is_root()) + { + ofs_running << " @IPI DRIVER: CHGDATA nelect=" << std::setprecision(17) + << charge[0] << " flag=" << flags << std::endl; + } +} + namespace { struct PosdataPayload @@ -177,18 +245,23 @@ struct PosdataPayload // 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 read_posdata(IpiSocket& socket, const DriverContext& context) { PosdataPayload payload; int io_failed = 0; std::string io_message; if (is_root()) { - if (state != DriverState::Ready) + if (context.state != DriverState::Ready) { io_failed = 1; io_message = "POSDATA requires READY state"; } + else if (context.cp.enabled && !context.charge_pending) + { + io_failed = 1; + io_message = "constant-potential POSDATA requires one fresh CHGDATA"; + } else { try @@ -212,7 +285,7 @@ PosdataPayload read_posdata(IpiSocket& socket, const UnitCell& ucell, const Driv std::size_t coordinate_count = 0; if (io_failed == 0 && !SocketFrame::checked_position_count(payload.nat_socket, - ucell.nat, + context.ucell->nat, coordinate_count, io_message)) { @@ -259,18 +332,29 @@ void bcast_posdata(PosdataPayload& payload) bcast_double_vector(payload.positions); } -void check_posdata_geometry(const DriverContext& context, const PosdataPayload& payload, DriverContext& mutable_context) +void check_posdata_geometry(DriverContext& context, const PosdataPayload& payload) { - 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)) + const std::vector wire_cell(payload.cell.begin(), payload.cell.end()); + const auto& expected_cell = context.cp.enabled && !context.cp_wire_cell.empty() + ? context.cp_wire_cell : context.reference_cell; + double cell_tolerance = unchanged_cell_tolerance(payload.cell); + if (context.cp.enabled && context.cp_wire_cell.empty()) + { + // Allow the two applications' rounded length constants once; lock later wire cells. + double scale = 1.0; + for (const double value : context.reference_cell) { scale = std::max(scale, std::abs(value)); } + cell_tolerance = 1.0e-6 * scale; + } + if (max_abs_delta(wire_cell, expected_cell) > cell_tolerance) { ModuleBase::WARNING_QUIT("ABACUS socket", "variable-cell socket updates are not supported yet."); } - if (!mutable_context.checked_initial_positions) + if (context.cp.enabled && context.cp_wire_cell.empty()) { context.cp_wire_cell = wire_cell; } + if (!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()) + 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; " @@ -280,14 +364,11 @@ void check_posdata_geometry(const DriverContext& context, const PosdataPayload& } } -void run_esolver_for_positions(UnitCell& ucell, - ModuleESolver::ESolver* esolver, - const std::vector& positions, - const int istep) +void run_esolver_for_positions(DriverContext& context, const std::vector& positions) { try { - set_positions_from_ipi_bohr(ucell, positions); + set_positions_from_ipi_bohr(*context.ucell, positions); } catch (const std::exception& exc) { @@ -299,7 +380,18 @@ void run_esolver_for_positions(UnitCell& ucell, } try { - esolver->runner(ucell, istep); + if (context.cp.enabled) + { + context.charge_pending = false; + const int functional_type = XC_Functional::get_func_type(); + SocketCP::run_electronic_step(*context.esolver, *context.ucell, *context.parameters, + *context.inp, *context.system, context.pending_control, + functional_type, context.istep); + } + else + { + context.esolver->runner(*context.ucell, context.istep); + } } catch (const std::exception& exc) { @@ -421,17 +513,21 @@ void compute_stress(UnitCell& ucell, ModuleESolver::ESolver* esolver, ComputedFr void handle_posdata(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running) { - PosdataPayload payload = read_posdata(socket, *context.ucell, context.state); + PosdataPayload payload = read_posdata(socket, context); bcast_posdata(payload); - check_posdata_geometry(context, payload, context); + check_posdata_geometry(context, payload); #ifdef __JSON // Each POSDATA frame is one electronic-structure output step. Json::init_output_array_obj(); #endif - run_esolver_for_positions(*context.ucell, context.esolver, payload.positions, context.istep); + run_esolver_for_positions(context, payload.positions); ComputedFrame computed; computed.scf_converged = all_ranks_converged(context.esolver->conv_esolver); + if (context.cp.enabled && !computed.scf_converged) + { + throw std::runtime_error("SCF did not converge; refusing to return constant-potential results"); + } if (!computed.scf_converged && is_root()) { ModuleBase::WARNING("ABACUS socket", @@ -446,13 +542,30 @@ void handle_posdata(IpiSocket& socket, DriverContext& context, std::ofstream& of { compute_stress(*context.ucell, context.esolver, computed); } + if (context.cp.enabled) + { + try + { + computed.electronic = SocketCP::collect_observables(context.cp, *context.esolver, + *context.ucell, *context.parameters); + } + catch (const std::exception& exc) + { + fail_during_collective_stage("electronic_observables", exc.what()); + } + catch (...) + { + fail_during_collective_stage("electronic_observables", "unknown electronic observables failure"); + } + computed.electronic_present = true; + } computed.valid = true; context.published = computed; ++context.istep; context.state = DriverState::HasData; } -void handle_getforce(IpiSocket& socket, DriverContext& context) +void handle_getforce(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running) { int io_failed = 0; std::string io_message; @@ -480,6 +593,7 @@ void handle_getforce(IpiSocket& socket, DriverContext& context) } socket.write_int32(static_cast(extra.size())); socket.write_string(extra); + ofs_running << " @IPI DRIVER: Sending extras JSON = " << extra << std::endl; } catch (const std::exception& exc) { @@ -492,8 +606,12 @@ void handle_getforce(IpiSocket& socket, DriverContext& context) context.state = DriverState::Ready; } -void handle_exit(std::ofstream& ofs_running) +void handle_exit(const DriverContext& context, std::ofstream& ofs_running) { + const bool unfinished = context.cp.enabled + && (context.state == DriverState::HasData || context.charge_pending); + const int local_failed = unfinished ? 1 : 0; + throw_if_any_rank_failed(local_failed, "EXIT with an unfinished request/result"); if (is_root()) { ofs_running << " ABACUS socket driver received i-PI EXIT" << std::endl; diff --git a/source/source_relax/socket_handlers.h b/source/source_relax/socket_handlers.h index f73224503ca..f69ab114985 100644 --- a/source/source_relax/socket_handlers.h +++ b/source/source_relax/socket_handlers.h @@ -9,6 +9,8 @@ #include class UnitCell; +class Parameter; +struct System_para; namespace ModuleESolver { @@ -23,28 +25,36 @@ struct DriverContext ModuleESolver::ESolver* esolver = nullptr; UnitCell* ucell = nullptr; const Input_para* inp = nullptr; + Parameter* parameters = nullptr; + const System_para* system = 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; + SocketCP::Configuration cp; + elecstate::ElectronicControl pending_control; + bool charge_pending = false; + std::vector cp_wire_cell; }; // 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); +std::string read_header_bcast(IpiSocket& socket, const DriverContext& context); void handle_status(IpiSocket& socket, const SocketUtils::DriverState state); void handle_init(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running); +void handle_chgdata(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_getforce(IpiSocket& socket, DriverContext& context, std::ofstream& ofs_running); -void handle_exit(std::ofstream& ofs_running); +void handle_exit(const DriverContext& context, std::ofstream& ofs_running); } // namespace SocketHandlers #endif diff --git a/source/source_relax/socket_utils.cpp b/source/source_relax/socket_utils.cpp index 8f41d005460..605941a135b 100644 --- a/source/source_relax/socket_utils.cpp +++ b/source/source_relax/socket_utils.cpp @@ -77,7 +77,14 @@ std::string properties_extra(const ComputedFrame& frame) { extra << ",\"stress\""; } - extra << "],\"scf_converged\":" << (frame.scf_converged ? "true" : "false") << "}"; + extra << "],\"scf_converged\":" + << (frame.scf_converged ? "true" : "false"); + if (frame.electronic_present) + { + const std::string electronic = SocketCP::serialize_observables(frame.electronic); + extra << "," << electronic.substr(1, electronic.size() - 2); + } + extra << "}"; return extra.str(); } @@ -288,6 +295,9 @@ void set_positions_from_ipi_bohr(UnitCell& ucell, const std::vector& pos 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->dis[ia].x -= std::round(atom->dis[ia].x); + atom->dis[ia].y -= std::round(atom->dis[ia].y); + atom->dis[ia].z -= std::round(atom->dis[ia].z); atom->taud[ia].x = dx; atom->taud[ia].y = dy; atom->taud[ia].z = dz; diff --git a/source/source_relax/socket_utils.h b/source/source_relax/socket_utils.h index 4314d165302..d69db6cc517 100644 --- a/source/source_relax/socket_utils.h +++ b/source/source_relax/socket_utils.h @@ -1,6 +1,7 @@ #ifndef ABACUS_SOURCE_RELAX_SOCKET_UTILS_H #define ABACUS_SOURCE_RELAX_SOCKET_UTILS_H +#include "source_relax/socket_constant_potential.h" #include "source_base/matrix.h" #include "source_relax/socket_frame.h" @@ -35,6 +36,8 @@ struct ComputedFrame bool forces_present = false; bool stress_present = false; bool scf_converged = true; + bool electronic_present = false; + SocketCP::Observables electronic; double energy_hartree = 0.0; std::vector forces_hartree_per_bohr; SocketFrame::Matrix9 virial_wire_hartree = {{0.0}}; diff --git a/source/source_relax/unittests/CMakeLists.txt b/source/source_relax/unittests/CMakeLists.txt index 2ab801f1dea..319ac4f5151 100644 --- a/source/source_relax/unittests/CMakeLists.txt +++ b/source/source_relax/unittests/CMakeLists.txt @@ -11,6 +11,10 @@ addtest(TARGET MODULE_RELAX_socket_frame_test SOURCES test_socket_frame.cpp ../socket_frame.cpp ../socket_frame_utils.cpp) +addtest(TARGET MODULE_RELAX_socket_constant_potential_test + SOURCES test_socket_constant_potential.cpp + ../socket_constant_potential.cpp) + set(socket_json_sources) if(ENABLE_JSON) list(APPEND socket_json_sources ../../source_io/module_json/abacusjson.cpp @@ -25,6 +29,7 @@ addtest( ../socket_driver.cpp ../socket_utils.cpp ../socket_handlers.cpp + ../socket_constant_potential.cpp ../socket_frame_utils.cpp ../socket_frame.cpp ../socket_ipi.cpp diff --git a/source/source_relax/unittests/test_relax_driver.cpp b/source/source_relax/unittests/test_relax_driver.cpp index 3f832eabf58..ba017103ec2 100644 --- a/source/source_relax/unittests/test_relax_driver.cpp +++ b/source/source_relax/unittests/test_relax_driver.cpp @@ -6,6 +6,7 @@ #include "source_relax/relax_driver.h" #include "source_relax/relax_stru_io.h" #include "source_relax/socket_driver.h" +#include "source_io/module_parameter/parameter.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -105,8 +106,22 @@ void write_stru(UnitCell&, } } // namespace relax_stru_io -void Socket_Driver::socket_driver(ModuleESolver::ESolver*, UnitCell&, const Input_para&, std::ofstream&) +class TestParameters { + public: + static Input_para& input(Parameter& parameters) { return parameters.input; } +}; + +namespace +{ +int socket_frame_count = 0; +} + +bool Socket_Driver::socket_driver(ModuleESolver::ESolver* solver, UnitCell& cell, + Parameter&, const Input_para&, const System_para&, std::ofstream&) +{ + for (int step = 0; step < socket_frame_count; ++step) { solver->runner(cell, step); } + return socket_frame_count > 0; } namespace @@ -179,6 +194,7 @@ class RelaxDriverTest : public ::testing::Test void SetUp() override { log_file_ = testing::TempDir() + "running.log"; + socket_frame_count = 0; } static void make_cell(UnitCell& ucell) @@ -212,7 +228,8 @@ class RelaxDriverTest : public ::testing::Test TEST_F(RelaxDriverTest, RelaxNmaxZeroDoesNotRunESolver) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 0; @@ -224,7 +241,7 @@ TEST_F(RelaxDriverTest, RelaxNmaxZeroDoesNotRunESolver) std::ofstream ofs(log_file_); Relax_Driver driver; testing::internal::CaptureStdout(); - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); const std::string output = testing::internal::GetCapturedStdout(); ofs.close(); @@ -235,7 +252,8 @@ TEST_F(RelaxDriverTest, RelaxNmaxZeroDoesNotRunESolver) TEST_F(RelaxDriverTest, FirstStepConvergesOnNonSimultaneousPath) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 3; @@ -250,7 +268,7 @@ TEST_F(RelaxDriverTest, FirstStepConvergesOnNonSimultaneousPath) std::ofstream ofs(log_file_); Relax_Driver driver; testing::internal::CaptureStdout(); - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); const std::string output = testing::internal::GetCapturedStdout(); ofs.close(); @@ -262,7 +280,8 @@ TEST_F(RelaxDriverTest, FirstStepConvergesOnNonSimultaneousPath) TEST_F(RelaxDriverTest, FirstStepConvergesOnSimultaneousPath) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "cg"; inp.relax_nmax = 3; @@ -277,7 +296,7 @@ TEST_F(RelaxDriverTest, FirstStepConvergesOnSimultaneousPath) std::ofstream ofs(log_file_); Relax_Driver driver; testing::internal::CaptureStdout(); - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); const std::string output = testing::internal::GetCapturedStdout(); ofs.close(); @@ -289,7 +308,8 @@ TEST_F(RelaxDriverTest, FirstStepConvergesOnSimultaneousPath) TEST_F(RelaxDriverTest, MaxStepsMovesFinalNonSimultaneousGeometry) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 1; @@ -304,7 +324,7 @@ TEST_F(RelaxDriverTest, MaxStepsMovesFinalNonSimultaneousGeometry) DriverTestESolver esolver(1.0); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); EXPECT_EQ(esolver.runner_count, 1); @@ -315,7 +335,8 @@ TEST_F(RelaxDriverTest, MaxStepsMovesFinalNonSimultaneousGeometry) TEST_F(RelaxDriverTest, MaxStepsStopsOnSimultaneousPath) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "cg"; inp.relax_nmax = 1; @@ -329,7 +350,7 @@ TEST_F(RelaxDriverTest, MaxStepsStopsOnSimultaneousPath) DriverTestESolver esolver(1.0); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); EXPECT_EQ(esolver.runner_count, 1); @@ -342,7 +363,8 @@ TEST_F(RelaxDriverTest, MaxStepsAfterAllUsesEvaluatedCellRelaxGeometry) { unitcell::allow_cell_change = true; - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "cell-relax"; inp.relax_method = "cg"; inp.relax_nmax = 1; @@ -361,7 +383,7 @@ TEST_F(RelaxDriverTest, MaxStepsAfterAllUsesEvaluatedCellRelaxGeometry) esolver.set_stress(1.0); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); unitcell::allow_cell_change = false; @@ -376,7 +398,8 @@ TEST_F(RelaxDriverTest, MaxStepsAfterAllUsesEvaluatedCellRelaxGeometry) TEST_F(RelaxDriverTest, FinalAllowedStepConvergesAndCallsAfterAllOnce) { - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 1; @@ -390,7 +413,7 @@ TEST_F(RelaxDriverTest, FinalAllowedStepConvergesAndCallsAfterAllOnce) DriverTestESolver esolver(1.0e-4); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); EXPECT_EQ(esolver.runner_count, 1); @@ -408,7 +431,8 @@ TEST_F(RelaxDriverTest, ExitStopsBeforeNonSimultaneousRelaxStep) exit_file << "stop_ion true\n"; } - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "lbfgs"; inp.relax_nmax = 3; @@ -422,7 +446,7 @@ TEST_F(RelaxDriverTest, ExitStopsBeforeNonSimultaneousRelaxStep) DriverTestESolver esolver(1.0); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); ASSERT_EQ(::chdir(old_cwd), 0); @@ -443,7 +467,8 @@ TEST_F(RelaxDriverTest, ExitStopsBeforeSimultaneousRelaxStep) exit_file << "stop_ion true\n"; } - Input_para inp; + Parameter parameters; + Input_para& inp = TestParameters::input(parameters); inp.calculation = "relax"; inp.relax_method = "cg"; inp.relax_nmax = 3; @@ -457,7 +482,7 @@ TEST_F(RelaxDriverTest, ExitStopsBeforeSimultaneousRelaxStep) DriverTestESolver esolver(1.0); std::ofstream ofs(log_file_); Relax_Driver driver; - driver.relax_driver(&esolver, ucell, inp, ofs); + driver.relax_driver(&esolver, ucell, parameters, ofs); ofs.close(); ASSERT_EQ(::chdir(old_cwd), 0); @@ -468,3 +493,32 @@ TEST_F(RelaxDriverTest, ExitStopsBeforeSimultaneousRelaxStep) } } // namespace + +TEST_F(RelaxDriverTest, SocketWithoutCompletedFrameSkipsFinalization) +{ + Parameter parameters; + TestParameters::input(parameters).socket_driver = true; + UnitCell cell; + make_cell(cell); + DriverTestESolver solver(0.0); + std::ofstream output(log_file_); + Relax_Driver driver; + driver.relax_driver(&solver, cell, parameters, output); + EXPECT_EQ(solver.runner_count, 0); + EXPECT_EQ(solver.after_all_count, 0); +} + +TEST_F(RelaxDriverTest, SocketCompletedFramesFinalizeOnce) +{ + Parameter parameters; + TestParameters::input(parameters).socket_driver = true; + socket_frame_count = 2; + UnitCell cell; + make_cell(cell); + DriverTestESolver solver(0.0); + std::ofstream output(log_file_); + Relax_Driver driver; + driver.relax_driver(&solver, cell, parameters, output); + EXPECT_EQ(solver.runner_count, 2); + EXPECT_EQ(solver.after_all_count, 1); +} diff --git a/source/source_relax/unittests/test_socket_constant_potential.cpp b/source/source_relax/unittests/test_socket_constant_potential.cpp new file mode 100644 index 00000000000..9d1e12a3a55 --- /dev/null +++ b/source/source_relax/unittests/test_socket_constant_potential.cpp @@ -0,0 +1,127 @@ +#include "source_relax/socket_constant_potential.h" + +#include "gtest/gtest.h" + +#include +#include +#include +#include + +namespace +{ +const std::string metadata = "NELECT:8.25,neutral_electrons:8"; +} + +TEST(SocketConstantPotentialTest, OrdinaryInitDoesNotEnableElectronControl) +{ + EXPECT_FALSE(SocketCP::parse_init("").enabled); + EXPECT_FALSE(SocketCP::parse_init("label:ordinary,unused:123").enabled); + EXPECT_THROW(SocketCP::make_control(SocketCP::Configuration(), 8.0, 1), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, FractionalElectronRequestIsPreserved) +{ + const auto config = SocketCP::parse_init(metadata); + const auto control = SocketCP::make_control(config, 8.3125, 0); + EXPECT_TRUE(config.enabled); + EXPECT_DOUBLE_EQ(8.3125, control.nelect); +} + +TEST(SocketConstantPotentialTest, ConstantPotentialInitRequiresCompleteMetadata) +{ + EXPECT_THROW(SocketCP::parse_init("NELECT:8"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("neutral_electrons:8"), std::runtime_error); + EXPECT_TRUE(SocketCP::parse_init("NELECT:8,neutral_electrons:8").enabled); + EXPECT_THROW(SocketCP::parse_init(metadata + ",NELECT:9"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init(metadata + ",neutral_electrons:8"), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, InitRejectsInvalidNumericalValues) +{ + EXPECT_THROW(SocketCP::parse_init("NELECT:8junk,neutral_electrons:8"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("NELECT:0,neutral_electrons:8"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("NELECT:8,neutral_electrons:-1"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("NELECT:8,neutral_electrons:8,solvation_update_stride:0"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("NELECT:8,neutral_electrons:8,solvation_update_stride:1.5"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init("NELECT:8,neutral_electrons:8,solvation_update_stride:2147483648"), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, PotentialWindowsSupportEachCartesianAxis) +{ + for (int axis = 0; axis < 3; ++axis) + { + const auto config = SocketCP::parse_init(metadata + ",potential_average_axis:" + + std::string(1, 'x' + axis) + + ",potential_average_min_A:12.5,potential_average_max_A:15"); + EXPECT_EQ(axis, config.window.axis); + EXPECT_DOUBLE_EQ(12.5, config.window.minimum_A); + EXPECT_DOUBLE_EQ(15.0, config.window.maximum_A); + } +} + +TEST(SocketConstantPotentialTest, PotentialWindowRejectsInvalidOrMissingBounds) +{ + EXPECT_THROW(SocketCP::parse_init(metadata + ",potential_average_axis:xy"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init(metadata + ",potential_average_axis:z,potential_average_min_A:1"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init(metadata + ",potential_average_axis:z,potential_average_min_A:2,potential_average_max_A:1"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init(metadata + ",potential_average_axis:z,potential_average_min_A:-1,potential_average_max_A:2"), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, RejectsRemovedStrideEvenForOrdinaryInit) +{ + EXPECT_THROW(SocketCP::parse_init("solvation_update_stride:1"), std::runtime_error); + EXPECT_THROW(SocketCP::parse_init(metadata + ",solvation_update_stride:5"), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, RejectsNeElectrodeAndUnknownFlags) +{ + const auto config = SocketCP::parse_init(metadata); + EXPECT_THROW(SocketCP::make_control(config, 8.0, 1), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, 8.0, 2), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, 8.0, 3), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, 8.0, 4), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, 8.0, -1), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, RejectsInvalidElectronRequests) +{ + const auto config = SocketCP::parse_init(metadata); + EXPECT_THROW(SocketCP::make_control(config, 0.0, 0), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, -1.0, 0), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, std::numeric_limits::quiet_NaN(), 0), std::runtime_error); + EXPECT_THROW(SocketCP::make_control(config, std::numeric_limits::infinity(), 0), std::runtime_error); +} + +TEST(SocketConstantPotentialTest, SerializesTheExistingElectronicExtrasKeys) +{ + SocketCP::Observables values; + values.nelect = 8.125; + values.fermi_eV = -0.25; + EXPECT_EQ("{\"fermi_level_eV\":-0.25,\"nelect\":8.125}", + SocketCP::serialize_observables(values)); + values.has_workfunction = true; + values.reference_eV = 4.0; + values.workfunction_eV = 4.25; + EXPECT_EQ("{\"fermi_level_eV\":-0.25,\"nelect\":8.125,\"workfunction_eV\":4.25," + "\"electrostatic_reference_eV\":4}", SocketCP::serialize_observables(values)); +} + +namespace +{ +class DecimalComma : public std::numpunct +{ + char do_decimal_point() const override { return ','; } +}; +} + +TEST(SocketConstantPotentialTest, JsonNumbersIgnoreTheGlobalDecimalSeparator) +{ + SocketCP::Observables values; + values.nelect = 8.125; + values.fermi_eV = -0.25; + const std::locale previous = std::locale(); + std::locale::global(std::locale(previous, new DecimalComma)); + const std::string json = SocketCP::serialize_observables(values); + std::locale::global(previous); + EXPECT_EQ("{\"fermi_level_eV\":-0.25,\"nelect\":8.125}", json); +} diff --git a/source/source_relax/unittests/test_socket_driver.cpp b/source/source_relax/unittests/test_socket_driver.cpp index 50e5fe25597..4e988642fe8 100644 --- a/source/source_relax/unittests/test_socket_driver.cpp +++ b/source/source_relax/unittests/test_socket_driver.cpp @@ -1,9 +1,17 @@ #include "for_test.h" +#ifdef __MPI #include "mpi.h" +#endif +#include "source_base/constants.h" #include "source_cell/unitcell.h" #include "source_esolver/esolver.h" +#include "source_hamilt/module_xc/xc_functional.h" #include "source_io/module_parameter/input_parameter.h" +#include "source_io/module_parameter/parameter.h" +#include "source_relax/socket_constant_potential.h" #include "source_relax/socket_driver.h" +#include "source_relax/socket_handlers.h" +#include "source_relax/socket_utils.h" #include "gmock/gmock.h" #include "gtest/gtest.h" @@ -24,6 +32,14 @@ #include #include +int XC_Functional::func_type = 2; + +class TestParameters +{ + public: + static Input_para& input(Parameter& parameters) { return parameters.input; } +}; + namespace { constexpr std::size_t kIpiHeaderLen = 12; @@ -244,6 +260,68 @@ class FakeESolver : public ModuleESolver::ESolver bool converged_; }; +class ElectronicObservablesSolver : public FakeESolver +{ + public: + ElectronicObservablesSolver() : FakeESolver(true) {} + + bool electronic_observables(const UnitCell&, const Parameter&, + const elecstate::PotentialWindow& window, + elecstate::ElectronicObservables& result) const override + { + ++collection_count; + last_window = window; + result = physical; + return supported; + } + + elecstate::ElectronicObservables physical; + bool supported = true; + mutable int collection_count = 0; + mutable elecstate::PotentialWindow last_window; +}; + +class ControlledESolver : public FakeESolver +{ + public: + ControlledESolver() : FakeESolver(true) {} + + bool apply_electronic_control(Parameter&, const elecstate::ElectronicControl& control) override + { + nelect = control.nelect; + ++control_count; + return true; + } + + void runner(BaseCell& cell, const int step) override + { + if (control_count != runner_count + 1) + { + throw std::runtime_error("each runner needs one electron control update"); + } + FakeESolver::runner(cell, step); + ++runner_count; + conv_esolver = true; + } + + double cal_energy() override { return FakeESolver::cal_energy() + nelect; } + + bool electronic_observables(const UnitCell&, const Parameter&, + const elecstate::PotentialWindow&, + elecstate::ElectronicObservables& result) const override + { + ++collection_count; + result.nelect = nelect; + result.fermi_ry = -0.125; + return true; + } + + double nelect = 0.0; + int control_count = 0; + int runner_count = 0; + mutable int collection_count = 0; +}; + struct DriverResult { int exit_code = -1; @@ -379,10 +457,8 @@ std::string read_pipe(const int fd) return output; } -DriverResult run_driver_frame(const bool converged, - const bool cal_force, - const bool cal_stress, - const std::function& peer_action) +DriverResult run_driver_session(const std::function& driver_action, + const std::function& peer_action) { UnixSocketServer server; int output_pipe[2]; @@ -406,15 +482,7 @@ DriverResult run_driver_frame(const bool converged, ::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); + driver_action(); std::cout.flush(); std::cerr.flush(); ::_exit(0); @@ -434,7 +502,6 @@ DriverResult run_driver_frame(const bool converged, { throw std::runtime_error(errno_message("setsockopt failed")); } - send_fixed_cell_frame(peer_fd); peer_action(peer_fd); } catch (...) @@ -467,8 +534,239 @@ DriverResult run_driver_frame(const bool converged, } return result; } + +DriverResult run_driver_frame(const bool converged, + const bool cal_force, + const bool cal_stress, + const std::function& peer_action) +{ + return run_driver_session([&]() { + UnitCell cell; + initialize_one_atom_cell(cell); + 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, cell, input, running); + }, [&](const int fd) { + send_fixed_cell_frame(fd); + peer_action(fd); + }); +} + +DriverResult run_constant_potential_session(const std::function& peer_action) +{ + return run_driver_session([]() { + UnitCell cell; + initialize_one_atom_cell(cell); + Parameter parameters; + Input_para& input = TestParameters::input(parameters); + input.cal_force = true; + input.cal_stress = false; + input.symmetry = "0"; + input.esolver_type = "ksdft"; + input.basis_type = "pw"; + input.device = "cpu"; + input.precision = "double"; + input.nspin = 1; + input.nbands = 8; + input.smearing_method = "gaussian"; + ControlledESolver solver; + std::ofstream running("/dev/null"); + Socket_Driver driver; + driver.socket_driver(&solver, cell, parameters, parameters.inp, parameters.globalv, running); + std::cout << "CP calls: " << solver.control_count << " controls, " + << solver.runner_count << " runners, " << solver.collection_count + << " collections\n"; + }, [&](const int fd) { + send_header(fd, "STATUS"); + EXPECT_EQ("NEEDINIT", read_header_or_close(fd)); + const std::string metadata = "NELECT:8,neutral_electrons:8"; + const std::int32_t replica = 0; + const std::int32_t nbytes = static_cast(metadata.size()); + send_header(fd, "INIT"); + send_value(fd, replica); + send_value(fd, nbytes); + send_all(fd, metadata.data(), metadata.size()); + peer_action(fd); + }); +} } // namespace +TEST(SocketDriverTest, RejectedElectronicStepLeavesSolverUntouched) +{ + UnitCell cell; + initialize_one_atom_cell(cell); + Parameter parameters; + // Configuration comes from the explicit references, not the opaque handle. + Input_para input; + System_para system; + input.esolver_type = "ksdft"; + input.basis_type = "pw"; + input.device = "cpu"; + input.precision = "double"; + input.nspin = 1; + input.nbands = 8; + input.smearing_method = "gaussian"; + ControlledESolver solver; + elecstate::ElectronicControl request; + request.nelect = 8.25; + + input.assume_isolated = "makov-payne"; + try + { + SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0); + FAIL() << "Makov-Payne request reached the solver"; + } + catch (const std::runtime_error& error) + { + EXPECT_THAT(error.what(), testing::HasSubstr("energy-only Makov-Payne")); + } + input.assume_isolated = "none"; + input.nbands = 4; + EXPECT_THROW(SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0), std::runtime_error); + input.nbands = 8; + input.smearing_method = "fixed"; + EXPECT_THROW(SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0), std::runtime_error); + input.smearing_method = "gaussian"; + system.use_uspp = true; + EXPECT_THROW(SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0), std::runtime_error); + system.use_uspp = false; + system.two_fermi = true; + EXPECT_THROW(SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0), std::runtime_error); + system.two_fermi = false; + EXPECT_THROW(SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 0, 0), std::runtime_error); + EXPECT_EQ(0, solver.control_count); + EXPECT_EQ(0, solver.runner_count); + EXPECT_DOUBLE_EQ(0.0, solver.nelect); + + input.smearing_method = "gaussian"; + SocketCP::run_electronic_step(solver, cell, parameters, input, system, request, 2, 0); + EXPECT_EQ(1, solver.control_count); + EXPECT_EQ(1, solver.runner_count); + EXPECT_DOUBLE_EQ(8.25, solver.nelect); + + FakeESolver unsupported(true); + EXPECT_THROW(SocketCP::run_electronic_step(unsupported, cell, parameters, input, system, request, 2, 0), std::runtime_error); +} + +TEST(SocketDriverTest, ConstantPotentialFramesApplyChargeBeforeEvaluation) +{ + const DriverResult result = run_constant_potential_session([](const int fd) { + const double counts[] = {8.25, 7.75}; + const double positions[] = {0.0, 0.125}; + const std::string extras[] = {"\"nelect\":8.25", "\"nelect\":7.75"}; + const std::int32_t flags = 0; + for (int frame = 0; frame < 2; ++frame) + { + send_header(fd, "STATUS"); + EXPECT_EQ("READY", read_header_or_close(fd)); + send_header(fd, "CHGDATA"); + send_value(fd, counts[frame]); + send_value(fd, flags); + send_header(fd, "STATUS"); + EXPECT_EQ("READY", read_header_or_close(fd)); + send_positions(fd, positions[frame]); + for (int query = 0; query < 2; ++query) + { + send_header(fd, "STATUS"); + EXPECT_EQ("HAVEDATA", read_header_or_close(fd)); + } + send_header(fd, "GETFORCE"); + const ForceResponse response = read_force_response(fd); + EXPECT_EQ("FORCEREADY", response.header); + EXPECT_DOUBLE_EQ((4.0 + positions[frame] + counts[frame]) * 0.5, response.energy_hartree); + ASSERT_EQ(3u, response.forces_hartree_per_bohr.size()); + EXPECT_DOUBLE_EQ(2.0 + positions[frame], response.forces_hartree_per_bohr[0]); + EXPECT_THAT(response.extra, testing::HasSubstr(extras[frame])); + EXPECT_THAT(response.extra, testing::HasSubstr("\"scf_converged\":true")); + } + send_header(fd, "STATUS"); + EXPECT_EQ("READY", read_header_or_close(fd)); + send_header(fd, "EXIT"); + EXPECT_EQ("", read_header_or_close(fd)); + }); + EXPECT_EQ(0, result.exit_code); + EXPECT_THAT(result.diagnostic, testing::HasSubstr("CP calls: 2 controls, 2 runners, 2 collections")); +} + +TEST(SocketDriverTest, ConstantPotentialPositionRequiresFreshCharge) +{ + const DriverResult result = run_constant_potential_session([](const int fd) { + send_header(fd, "POSDATA"); + EXPECT_EQ("", read_header_or_close(fd)); + }); + EXPECT_NE(0, result.exit_code); + EXPECT_THAT(result.diagnostic, testing::HasSubstr("constant-potential POSDATA requires one fresh CHGDATA")); +} + +TEST(SocketElectronicFrameTest, CollectionConvertsPhysicalUnitsAndKeepsZeroFermi) +{ + UnitCell cell; + Parameter parameters; + ElectronicObservablesSolver solver; + solver.physical.nelect = 8.125; + solver.physical.fermi_ry = 0.0; + solver.physical.reference_ry = 0.25; + const auto config = SocketCP::parse_init("NELECT:8,neutral_electrons:8," + "potential_average_axis:y,potential_average_min_A:2,potential_average_max_A:3"); + const auto values = SocketCP::collect_observables(config, solver, cell, parameters); + EXPECT_DOUBLE_EQ(0.0, values.fermi_eV); + EXPECT_DOUBLE_EQ(8.125, values.nelect); + EXPECT_TRUE(values.has_workfunction); + EXPECT_DOUBLE_EQ(0.25 * ModuleBase::Ry_to_eV, values.reference_eV); + EXPECT_DOUBLE_EQ(values.reference_eV, values.workfunction_eV); + EXPECT_EQ(1, solver.last_window.axis); + EXPECT_DOUBLE_EQ(2.0, solver.last_window.minimum_A); + EXPECT_DOUBLE_EQ(3.0, solver.last_window.maximum_A); + EXPECT_EQ(1, solver.collection_count); +} + +TEST(SocketElectronicFrameTest, CollectionRejectsMissingOrInvalidPhysicalObservables) +{ + UnitCell cell; + Parameter parameters; + ElectronicObservablesSolver solver; + const auto config = SocketCP::parse_init("NELECT:8,neutral_electrons:8"); + EXPECT_THROW(SocketCP::collect_observables(SocketCP::Configuration(), solver, cell, parameters), + std::runtime_error); + EXPECT_EQ(0, solver.collection_count); + solver.supported = false; + EXPECT_THROW(SocketCP::collect_observables(config, solver, cell, parameters), std::runtime_error); + solver.supported = true; + solver.physical.nelect = 0.0; + EXPECT_THROW(SocketCP::collect_observables(config, solver, cell, parameters), std::runtime_error); + solver.physical.nelect = 8.0; + solver.physical.fermi_ry = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(SocketCP::collect_observables(config, solver, cell, parameters), std::runtime_error); +} + +TEST(SocketElectronicFrameTest, SerializationUsesPublishedSnapshotWithoutCollectingAgain) +{ + UnitCell cell; + Parameter parameters; + ElectronicObservablesSolver solver; + solver.physical.nelect = 8.125; + solver.physical.fermi_ry = -0.125; + const auto config = SocketCP::parse_init("NELECT:8,neutral_electrons:8"); + SocketUtils::ComputedFrame frame; + frame.valid = true; + frame.forces_present = true; + frame.electronic = SocketCP::collect_observables(config, solver, cell, parameters); + frame.electronic_present = true; + solver.physical.nelect = 9.0; + solver.supported = false; + const std::string first = SocketUtils::properties_extra(frame); + const std::string second = SocketUtils::properties_extra(frame); + EXPECT_EQ(first, second); + EXPECT_THAT(first, testing::HasSubstr("\"nelect\":8.125")); + EXPECT_THAT(first, testing::HasSubstr("\"present\":[\"energy\",\"forces\"]")); + EXPECT_THAT(first, testing::Not(testing::HasSubstr("workfunction_eV"))); + EXPECT_EQ(1, solver.collection_count); +} + TEST(SocketDriverTest, NonconvergedFrameIsPublishedWithMetadata) { ForceResponse response; @@ -579,6 +877,31 @@ TEST(SocketDriverTest, ConsumedFrameCannotBeReturnedTwice) EXPECT_THAT(result.diagnostic, testing::HasSubstr("GETFORCE requires HAVEDATA")); } +TEST(SocketDriverTest, OrdinaryExitMayDiscardUnconsumedFrame) +{ + const DriverResult result = run_driver_frame(true, true, false, [&](const int fd) { + send_header(fd, "STATUS"); + EXPECT_EQ("HAVEDATA", read_header_or_close(fd)); + send_header(fd, "EXIT"); + EXPECT_EQ("", read_header_or_close(fd)); + }); + EXPECT_EQ(0, result.exit_code); +} + +TEST(SocketDriverTest, ConstantPotentialExitRequiresCompletedRequest) +{ + SocketHandlers::DriverContext context; + std::ofstream running("/dev/null"); + context.cp.enabled = true; + context.state = SocketUtils::DriverState::HasData; + EXPECT_THROW(SocketHandlers::handle_exit(context, running), std::runtime_error); + context.state = SocketUtils::DriverState::Ready; + context.charge_pending = true; + EXPECT_THROW(SocketHandlers::handle_exit(context, running), std::runtime_error); + context.charge_pending = false; + EXPECT_NO_THROW(SocketHandlers::handle_exit(context, running)); +} + TEST(SocketDriverTest, InvalidNextFrameCannotReturnPreviousResults) { const DriverResult result = run_driver_frame(true, true, false, [&](const int fd) { @@ -593,9 +916,13 @@ TEST(SocketDriverTest, InvalidNextFrameCannotReturnPreviousResults) int main(int argc, char** argv) { +#ifdef __MPI MPI_Init(&argc, &argv); +#endif testing::InitGoogleTest(&argc, argv); const int result = RUN_ALL_TESTS(); +#ifdef __MPI MPI_Finalize(); +#endif return result; }