From f5043e62e396a9c1071c1686673d1407b8890cff Mon Sep 17 00:00:00 2001 From: Kai Luo Date: Wed, 13 May 2026 11:00:06 +0800 Subject: [PATCH 1/3] Rename Makov-Payne energy term --- source/Makefile.Objects | 1 + source/source_estate/CMakeLists.txt | 1 + source/source_estate/elecstate_energy.cpp | 30 ++ source/source_estate/fp_energy.cpp | 7 +- source/source_estate/fp_energy.h | 1 + source/source_estate/makov_payne.cpp | 346 ++++++++++++++++++ source/source_estate/makov_payne.h | 29 ++ .../source_estate/module_pot/potential_new.h | 4 + .../module_parameter/input_parameter.h | 1 + .../read_input_item_system.cpp | 36 ++ 10 files changed, 453 insertions(+), 3 deletions(-) create mode 100644 source/source_estate/makov_payne.cpp create mode 100644 source/source_estate/makov_payne.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 803cd04018f..c7bda34fa96 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -232,6 +232,7 @@ OBJS_ELECSTAT=elecstate.o\ elecstate_pw.o\ elecstate_pw_sdft.o\ elecstate_pw_cal_tau.o\ + makov_payne.o\ elecstate_op.o\ efield.o\ gatefield.o\ diff --git a/source/source_estate/CMakeLists.txt b/source/source_estate/CMakeLists.txt index b1d4d9014f0..3391c5acf38 100644 --- a/source/source_estate/CMakeLists.txt +++ b/source/source_estate/CMakeLists.txt @@ -8,6 +8,7 @@ list(APPEND objects elecstate_pw.cpp elecstate_pw_sdft.cpp elecstate_pw_cal_tau.cpp + makov_payne.cpp module_pot/gatefield.cpp module_pot/efield.cpp module_pot/H_Hartree_pw.cpp diff --git a/source/source_estate/elecstate_energy.cpp b/source/source_estate/elecstate_energy.cpp index 3c58360b902..34e6c3aa278 100644 --- a/source/source_estate/elecstate_energy.cpp +++ b/source/source_estate/elecstate_energy.cpp @@ -2,7 +2,9 @@ #include "source_base/global_variable.h" #include "source_base/parallel_comm.h" #include "source_base/parallel_reduce.h" +#include "makov_payne.h" #include "source_hamilt/module_xc/xc_functional.h" +#include "source_estate/module_pot/H_Hartree_pw.h" #include "source_io/module_parameter/parameter.h" #include @@ -340,6 +342,34 @@ void ElecState::cal_energies(const int type) this->f_en.e_local_pp = get_local_pp_energy(); + if (PARAM.inp.assume_isolated == "makov-payne") + { + const UnitCell* ucell = this->pot->get_ucell(); + if (ucell == nullptr || this->charge == nullptr || this->charge->rhopw == nullptr) + { + ModuleBase::WARNING_QUIT("ElecState::cal_energies", + "Makov-Payne correction requires an initialized unit cell and charge density."); + } + std::vector v_elecstat; + const double* v_elecstat_ptr = nullptr; + { + ModuleBase::matrix vh(PARAM.inp.nspin, this->charge->rhopw->nrxx); + vh = elecstate::H_Hartree_pw::v_hartree(*ucell, this->charge->rhopw, PARAM.inp.nspin, this->charge->rho); + v_elecstat.assign(this->charge->rhopw->nrxx, 0.0); + const double* v_fixed = this->pot->get_fixed_v(); + for (int ir = 0; ir < this->charge->rhopw->nrxx; ++ir) + { + v_elecstat[ir] = vh(0, ir) + v_fixed[ir]; + } + v_elecstat_ptr = v_elecstat.data(); + } + this->f_en.correction_el = makov_payne_correction(*ucell, *this->charge, v_elecstat_ptr).total; + } + else + { + this->f_en.correction_el = 0.0; + } + #ifdef __MLALGO this->f_en.ml_exx = this->pot->get_ml_exx_energy(); #endif diff --git a/source/source_estate/fp_energy.cpp b/source/source_estate/fp_energy.cpp index a8d3124b4f7..99236a634e5 100644 --- a/source/source_estate/fp_energy.cpp +++ b/source/source_estate/fp_energy.cpp @@ -16,7 +16,7 @@ namespace elecstate double fenergy::calculate_etot() { etot = eband + deband + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield - + gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf + escon + ml_exx; + + gatefield + evdw + correction_el + esol_el + esol_cav + edftu + edeepks_scf + escon + ml_exx; return etot; } @@ -24,7 +24,7 @@ double fenergy::calculate_etot() double fenergy::calculate_harris() { etot_harris = eband + deband_harris + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx - + efield + gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf + escon + ml_exx; + + efield + gatefield + evdw + correction_el + esol_el + esol_cav + edftu + edeepks_scf + escon + ml_exx; return etot_harris; } @@ -32,7 +32,7 @@ double fenergy::calculate_harris() void fenergy::clear_all() { etot = etot_old = eband = deband = etxc = etxcc = vtxc = ewald_energy = hartree_energy = demet = descf = exx - = efield = gatefield = evdw = etot_harris = deband_harris = esol_el = esol_cav = edftu = edeepks_scf = escon + = efield = gatefield = evdw = correction_el = etot_harris = deband_harris = esol_el = esol_cav = edftu = edeepks_scf = escon = ml_exx = 0.0; } @@ -53,6 +53,7 @@ void fenergy::print_all() const std::cout << " efiled=" << efield << std::endl; std::cout << " gatefiled=" << gatefield << std::endl; std::cout << " evdw=" << evdw << std::endl; + std::cout << " correction_el=" << correction_el << std::endl; std::cout << " esol_el=" << esol_el << std::endl; std::cout << " esol_cav=" << esol_cav << std::endl; std::cout << " edftu=" << edftu << std::endl; diff --git a/source/source_estate/fp_energy.h b/source/source_estate/fp_energy.h index 3f37ca26563..5577947f44a 100644 --- a/source/source_estate/fp_energy.h +++ b/source/source_estate/fp_energy.h @@ -34,6 +34,7 @@ struct fenergy double efield = 0.0; ///< dipole potential in surface calculations double gatefield = 0.0; ///< correction energy for gatefield double evdw = 0.0; ///< the vdw energy + double correction_el = 0.0; ///< electrostatic isolated-cell correction double etot_harris = 0.0; ///< total energy of harris functional double deband_harris = 0.0; ///< correction for harris energy diff --git a/source/source_estate/makov_payne.cpp b/source/source_estate/makov_payne.cpp new file mode 100644 index 00000000000..7f18d3800c3 --- /dev/null +++ b/source/source_estate/makov_payne.cpp @@ -0,0 +1,346 @@ +#include "makov_payne.h" + +#include "source_base/constants.h" +#include "source_base/global_variable.h" +#include "source_base/parallel_reduce.h" +#include "source_base/tool_quit.h" +#include "source_cell/unitcell.h" +#include "source_estate/module_charge/charge.h" +#include "source_io/module_parameter/parameter.h" + +#include +#include +#include +#include +#include +#include + +namespace +{ +constexpr double madelung_sc = 2.8373; +constexpr double madelung_fcc = 2.8883; +constexpr double madelung_bcc = 2.8885; + +int cubic_ibrav(const std::string& latname) +{ + if (latname == "sc") + { + return 1; + } + if (latname == "fcc") + { + return 2; + } + if (latname == "bcc") + { + return 3; + } + return 0; +} + +double madelung_constant(const int ibrav) +{ + if (ibrav == 1) + { + return madelung_sc; + } + if (ibrav == 2) + { + return madelung_fcc; + } + return madelung_bcc; +} + +ModuleBase::Vector3 frac_to_cart_lat0(const UnitCell& ucell, const ModuleBase::Vector3& frac) +{ + return frac.x * ucell.a1 + frac.y * ucell.a2 + frac.z * ucell.a3; +} + +double wrap_delta(double x) +{ + x -= std::floor(x + 0.5); + return x; +} + +void add_average(const Charge& charge, const double* values, const int direction, const int length, std::vector& ave) +{ + const ModulePW::PW_Basis* rhopw = charge.rhopw; + for (int ir = 0; ir < rhopw->nrxx; ++ir) + { + int index = 0; + if (direction == 0) + { + index = ir / (rhopw->ny * rhopw->nplane); + } + else if (direction == 1) + { + const int ix = ir / (rhopw->ny * rhopw->nplane); + index = ir / rhopw->nplane - ix * rhopw->ny; + } + else + { + index = ir % rhopw->nplane + rhopw->startz_current; + } + ave[index] += values[ir]; + } +#ifdef __MPI + Parallel_Reduce::reduce_pool(ave.data(), length); +#endif + const int surface = rhopw->nxyz / length; + for (double& v : ave) + { + v /= static_cast(surface); + } +} + +bool estimate_corrected_vacuum_level(const UnitCell& ucell, + const Charge& charge, + const double* v_elecstat, + const ModuleBase::Vector3& x0_frac, + const double net_charge, + double& vacuum_level_ev) +{ + if (v_elecstat == nullptr || charge.rhopw == nullptr || charge.rhopw->nxyz <= 0) + { + return false; + } + + double vacuum[3] = {0.0, 0.0, 0.0}; + for (int dir = 0; dir < 3; ++dir) + { + std::vector pos; + pos.reserve(ucell.nat); + for (int it = 0; it < ucell.ntype; ++it) + { + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + pos.push_back(ucell.atoms[it].taud[ia][dir]); + } + } + std::sort(pos.begin(), pos.end()); + for (std::size_t i = 1; i < pos.size(); ++i) + { + vacuum[dir] = std::max(vacuum[dir], pos[i] - pos[i - 1]); + } + if (!pos.empty()) + { + vacuum[dir] = std::max(vacuum[dir], pos.front() + 1.0 - pos.back()); + } + } + + int direction = 2; + const double lengths[3] = {ucell.a1.norm() * ucell.lat0, ucell.a2.norm() * ucell.lat0, ucell.a3.norm() * ucell.lat0}; + vacuum[0] *= lengths[0]; + vacuum[1] *= lengths[1]; + vacuum[2] *= lengths[2]; + if (vacuum[0] > vacuum[direction]) + { + direction = 0; + } + if (vacuum[1] > vacuum[direction]) + { + direction = 1; + } + + const ModulePW::PW_Basis* rhopw = charge.rhopw; + int length = rhopw->nz; + if (direction == 0) + { + length = rhopw->nx; + } + else if (direction == 1) + { + length = rhopw->ny; + } + + std::vector total_rho(rhopw->nrxx, 0.0); + for (int ir = 0; ir < rhopw->nrxx; ++ir) + { + total_rho[ir] = std::fabs(charge.rho[0][ir]); + } + if (PARAM.inp.nspin == 2) + { + for (int ir = 0; ir < rhopw->nrxx; ++ir) + { + total_rho[ir] += std::fabs(charge.rho[1][ir]); + } + } + + std::vector rho_ave(length, 0.0); + add_average(charge, total_rho.data(), direction, length, rho_ave); + + int min_index = 0; + double min_value = 1.0e100; + const double windows[7] = {0.1, 0.2, 0.3, 0.4, 0.3, 0.2, 0.1}; + for (int i = 0; i < length; ++i) + { + double sum = 0.0; + const int temp = i - 3 + length; + for (int win = 0; win < 7; ++win) + { + sum += rho_ave[(temp + win) % length] * windows[win]; + } + if (sum < min_value) + { + min_value = sum; + min_index = i; + } + } + + std::vector v_ave(length, 0.0); + add_average(charge, v_elecstat, direction, length, v_ave); + + ModuleBase::Vector3 vac_frac = x0_frac; + vac_frac[direction] = (static_cast(min_index) + 0.5) / static_cast(length); + ModuleBase::Vector3 dfrac(wrap_delta(vac_frac.x - x0_frac.x), + wrap_delta(vac_frac.y - x0_frac.y), + wrap_delta(vac_frac.z - x0_frac.z)); + double r = frac_to_cart_lat0(ucell, dfrac).norm() * ucell.lat0; + if (r < 1.0e-12) + { + r = 0.5 * std::min(lengths[0], std::min(lengths[1], lengths[2])); + } + + vacuum_level_ev = (v_ave[min_index] + ModuleBase::e2 * net_charge / r) * ModuleBase::Ry_to_eV; + return true; +} +} // namespace + +namespace elecstate +{ + +MakovPayneResult makov_payne_correction(const UnitCell& ucell, + const Charge& charge, + const double* v_elecstat, + bool print) +{ + const int ibrav = cubic_ibrav(PARAM.inp.latname); + if (ibrav == 0) + { + ModuleBase::WARNING_QUIT("Makov-Payne", "Makov-Payne correction is defined only for cubic lattices: latname = sc, fcc, or bcc."); + } + if (charge.rhopw == nullptr || charge.rho == nullptr) + { + ModuleBase::WARNING_QUIT("Makov-Payne", "charge density is not available."); + } + + double zion = 0.0; + ModuleBase::Vector3 x0_frac(0.0, 0.0, 0.0); + for (int it = 0; it < ucell.ntype; ++it) + { + const double zv = ucell.atoms[it].ncpp.zv; + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + zion += zv; + x0_frac += ucell.atoms[it].taud[ia] * zv; + } + } + x0_frac = x0_frac / zion; + + ModuleBase::Vector3 dipole_ion(0.0, 0.0, 0.0); + ModuleBase::Vector3 quadrupole_ion(0.0, 0.0, 0.0); + for (int it = 0; it < ucell.ntype; ++it) + { + const double zv = ucell.atoms[it].ncpp.zv; + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) + { + ModuleBase::Vector3 dfrac(wrap_delta(ucell.atoms[it].taud[ia].x - x0_frac.x), + wrap_delta(ucell.atoms[it].taud[ia].y - x0_frac.y), + wrap_delta(ucell.atoms[it].taud[ia].z - x0_frac.z)); + ModuleBase::Vector3 dr = frac_to_cart_lat0(ucell, dfrac) * ucell.lat0; + dipole_ion += dr * zv; + quadrupole_ion.x += zv * dr.x * dr.x; + quadrupole_ion.y += zv * dr.y * dr.y; + quadrupole_ion.z += zv * dr.z * dr.z; + } + } + + const ModulePW::PW_Basis* rhopw = charge.rhopw; + const double dv = rhopw->omega / static_cast(rhopw->nxyz); + double electron_number = 0.0; + ModuleBase::Vector3 dipole_el(0.0, 0.0, 0.0); + ModuleBase::Vector3 quadrupole_el(0.0, 0.0, 0.0); + + for (int ir = 0; ir < rhopw->nrxx; ++ir) + { + const int ix = ir / (rhopw->ny * rhopw->nplane); + const int iy = ir / rhopw->nplane - ix * rhopw->ny; + const int iz = ir % rhopw->nplane + rhopw->startz_current; + double rho = charge.rho[0][ir]; + if (PARAM.inp.nspin == 2) + { + rho += charge.rho[1][ir]; + } + const double weight = rho * dv; + electron_number += weight; + ModuleBase::Vector3 frac((static_cast(ix) + 0.5) / static_cast(rhopw->nx), + (static_cast(iy) + 0.5) / static_cast(rhopw->ny), + (static_cast(iz) + 0.5) / static_cast(rhopw->nz)); + ModuleBase::Vector3 dfrac(wrap_delta(frac.x - x0_frac.x), + wrap_delta(frac.y - x0_frac.y), + wrap_delta(frac.z - x0_frac.z)); + ModuleBase::Vector3 dr = frac_to_cart_lat0(ucell, dfrac) * ucell.lat0; + dipole_el += dr * weight; + quadrupole_el.x += weight * dr.x * dr.x; + quadrupole_el.y += weight * dr.y * dr.y; + quadrupole_el.z += weight * dr.z * dr.z; + } + +#ifdef __MPI + Parallel_Reduce::reduce_pool(electron_number); + Parallel_Reduce::reduce_pool(&dipole_el.x, 3); + Parallel_Reduce::reduce_pool(&quadrupole_el.x, 3); +#endif + + MakovPayneResult result; + result.charge = zion - electron_number; + const ModuleBase::Vector3 dipole = dipole_ion - dipole_el; + const ModuleBase::Vector3 quadrupole = quadrupole_ion - quadrupole_el; + const double aa = quadrupole.x + quadrupole.y + quadrupole.z; + const double bb = dipole.x * dipole.x + dipole.y * dipole.y + dipole.z * dipole.z; + const double alpha = madelung_constant(ibrav); + + const double corr1 = -alpha / ucell.lat0 * result.charge * result.charge / 2.0 * ModuleBase::e2; + const double corr2 = (2.0 / 3.0 * ModuleBase::PI) * (result.charge * aa - bb) / std::pow(ucell.lat0, 3) + * ModuleBase::e2; + result.first_order = -corr1; + result.second_order = -corr2; + result.total = result.first_order + result.second_order; + result.has_vacuum_level = estimate_corrected_vacuum_level(ucell, + charge, + v_elecstat, + x0_frac, + result.charge, + result.vacuum_level); + + if (print && GlobalV::MY_RANK == 0) + { + GlobalV::ofs_running << std::setprecision(8) << std::fixed; + GlobalV::ofs_running << "\n ********* MAKOV-PAYNE CORRECTION *********" << std::endl; + GlobalV::ofs_running << " Makov-Payne correction with Madelung constant = " << std::setw(8) << alpha + << std::endl; + GlobalV::ofs_running << " Makov-Payne correction " << std::setw(14) << result.first_order + << " Ry = " << std::setw(10) << result.first_order * ModuleBase::Ry_to_eV + << " eV (1st order, 1/a0)" << std::endl; + GlobalV::ofs_running << " " << std::setw(14) << result.second_order + << " Ry = " << std::setw(10) << result.second_order * ModuleBase::Ry_to_eV + << " eV (2nd order, 1/a0^3)" << std::endl; + GlobalV::ofs_running << " " << std::setw(14) << result.total + << " Ry = " << std::setw(10) << result.total * ModuleBase::Ry_to_eV + << " eV (total)" << std::endl; + if (result.has_vacuum_level) + { + GlobalV::ofs_running << " Corrected vacuum level = " << std::setw(16) << result.vacuum_level << " eV" + << std::endl; + } + else + { + GlobalV::ofs_running << " Corrected vacuum level was not evaluated because the electrostatic potential is unavailable." + << std::endl; + } + } + + return result; +} + +} // namespace elecstate diff --git a/source/source_estate/makov_payne.h b/source/source_estate/makov_payne.h new file mode 100644 index 00000000000..80a1307ff2c --- /dev/null +++ b/source/source_estate/makov_payne.h @@ -0,0 +1,29 @@ +#ifndef MAKOV_PAYNE_H +#define MAKOV_PAYNE_H + +class UnitCell; +class Charge; + +namespace elecstate +{ + +struct MakovPayneResult +{ + double charge = 0.0; + double first_order = 0.0; + double second_order = 0.0; + double total = 0.0; + double vacuum_level = 0.0; + bool has_vacuum_level = false; +}; + +/// Calculate and print the Makov-Payne isolated-cell correction for cubic cells. +/// All energies are in Ry; the optional corrected vacuum level is in eV. +MakovPayneResult makov_payne_correction(const UnitCell& ucell, + const Charge& charge, + const double* v_elecstat = nullptr, + bool print = true); + +} // namespace elecstate + +#endif // MAKOV_PAYNE_H diff --git a/source/source_estate/module_pot/potential_new.h b/source/source_estate/module_pot/potential_new.h index ca79ca02f63..804674cd39f 100644 --- a/source/source_estate/module_pot/potential_new.h +++ b/source/source_estate/module_pot/potential_new.h @@ -175,6 +175,10 @@ class Potential : public PotBase { return this->rho_basis_; } + const UnitCell* get_ucell() const + { + return this->ucell_; + } // What about adding a function to get the wfc? // This is useful for the calculation of the exx energy diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 029ad364eb5..030d7c278f0 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -33,6 +33,7 @@ struct Input_para int kpar = 1; ///< ecch pool is for one k point int bndpar = 1; ///< parallel for stochastic/deterministic bands std::string latname = "user_defined_lattice"; ///< lattice name + std::string assume_isolated = "none"; ///< isolated-system correction: none or makov-payne double ecutwfc = 0; ///< energy cutoff for wavefunctions double ecutrho = 0; ///< energy cutoff for charge/potential diff --git a/source/source_io/module_parameter/read_input_item_system.cpp b/source/source_io/module_parameter/read_input_item_system.cpp index d746d2b5b12..651788b5080 100644 --- a/source/source_io/module_parameter/read_input_item_system.cpp +++ b/source/source_io/module_parameter/read_input_item_system.cpp @@ -384,6 +384,42 @@ Available options are: read_sync_string(input.latname); this->add_item(item); } + { + Input_Item item("assume_isolated"); + item.annotation = "isolated-system electrostatic correction"; + item.category = "System variables"; + item.type = "String"; + item.description = R"(Used to perform a calculation assuming an isolated system in a 3D supercell. + +Available options are: +* none: regular periodic calculation without isolated-system correction. +* makov-payne, m-p, mp: compute the Makov-Payne correction to the total energy and estimate a corrected vacuum level for eigenvalue alignment. This option is available only for cubic lattices (latname = sc, fcc, or bcc). + +Theory: G. Makov and M. C. Payne, Phys. Rev. B 51, 4014 (1995).)"; + item.default_value = "none"; + read_sync_string(input.assume_isolated); + item.reset_value = [](const Input_Item& item, Parameter& para) { + if (para.input.assume_isolated == "m-p" || para.input.assume_isolated == "mp") + { + para.input.assume_isolated = "makov-payne"; + } + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + const std::vector allowed = {"none", "makov-payne", "m-p", "mp"}; + if (std::find(allowed.begin(), allowed.end(), para.input.assume_isolated) == allowed.end()) + { + ModuleBase::WARNING_QUIT("ReadInput", nofound_str(allowed, "assume_isolated")); + } + if ((para.input.assume_isolated == "makov-payne" || para.input.assume_isolated == "m-p" + || para.input.assume_isolated == "mp") + && para.input.latname != "sc" && para.input.latname != "fcc" && para.input.latname != "bcc") + { + ModuleBase::WARNING_QUIT("ReadInput", + "Makov-Payne correction is available only for cubic lattices: latname = sc, fcc, or bcc."); + } + }; + this->add_item(item); + } { Input_Item item("init_wfc"); item.annotation = "start wave functions are from 'atomic', " From 827863194ea3eed507cbf979cf751996b9c8e329 Mon Sep 17 00:00:00 2001 From: Kai Luo Date: Mon, 18 May 2026 10:05:02 +0800 Subject: [PATCH 2/3] Fix elecstate energy unit test link dependencies --- source/source_estate/test/CMakeLists.txt | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index aa69f4e29d8..f8060373ed0 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -64,8 +64,12 @@ AddTest( AddTest( TARGET MODULE_ESTATE_elecstate_energy - LIBS parameter ${math_libs} base device - SOURCES elecstate_energy_test.cpp ../elecstate_energy.cpp ../fp_energy.cpp + LIBS parameter ${math_libs} base device planewave_serial + SOURCES elecstate_energy_test.cpp + ../elecstate_energy.cpp + ../fp_energy.cpp + ../makov_payne.cpp + ../module_pot/H_Hartree_pw.cpp ) AddTest( From 97a552b790ee80149812cb9875b4a6462d11babe Mon Sep 17 00:00:00 2001 From: Kai Luo Date: Tue, 19 May 2026 21:45:39 +0800 Subject: [PATCH 3/3] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- source/source_estate/makov_payne.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/source/source_estate/makov_payne.cpp b/source/source_estate/makov_payne.cpp index 7f18d3800c3..5796a954cc1 100644 --- a/source/source_estate/makov_payne.cpp +++ b/source/source_estate/makov_payne.cpp @@ -235,6 +235,11 @@ MakovPayneResult makov_payne_correction(const UnitCell& ucell, x0_frac += ucell.atoms[it].taud[ia] * zv; } } + if (zion <= 0.0) + { + ModuleBase::WARNING_QUIT("Makov-Payne", + "total ionic valence must be positive to compute the Makov-Payne correction."); + } x0_frac = x0_frac / zion; ModuleBase::Vector3 dipole_ion(0.0, 0.0, 0.0);