Electrostatic correction of Makov-Payne type

[kluonju]

With the help of codex, I have included the correction of Markov-Payne type for non-neutral isolated system. One can activate this function with assume_isolated mp and by default it is set none.

Note: this will only work with cubic system so far and we need to specify the lattice type using something like latname sc . This has been tested on LiH with fractional charge.

Further implementation on Martyna-Tuckerman type would be nice to have, which does not have such limitation.

[Copilot]

Pull request overview

Adds an isolated-supercell electrostatic correction option (Makov–Payne) for non-neutral systems, wired into the energy evaluation path and exposed via a new assume_isolated input keyword.

Changes:

• Introduces assume_isolated input parameter (default none) with Makov–Payne options and cubic-lattice validation.
• Implements Makov–Payne correction (+ optional corrected vacuum-level estimate) and integrates it into ElecState total-energy accounting via a new fenergy::correction_el term.
• Updates build/test targets and adds a Potential accessor to retrieve the associated UnitCell.

Reviewed changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
source/source_io/module_parameter/read_input_item_system.cpp	Adds assume_isolated input item, documentation, normalization, and validation.
source/source_io/module_parameter/input_parameter.h	Adds Input_para::assume_isolated with default none.
source/source_estate/elecstate_energy.cpp	Computes and applies Makov–Payne correction during energy evaluation.
source/source_estate/makov_payne.h	Declares Makov–Payne correction API and result struct.
source/source_estate/makov_payne.cpp	Implements Makov–Payne correction and vacuum-level estimation.
source/source_estate/fp_energy.h	Adds correction_el energy component.
source/source_estate/fp_energy.cpp	Includes correction_el in totals/clearing/printing.
source/source_estate/module_pot/potential_new.h	Adds Potential::get_ucell() accessor.
source/source_estate/CMakeLists.txt	Adds makov_payne.cpp to estate objects.
source/Makefile.Objects	Adds makov_payne.o to estate object list.
source/source_estate/test/CMakeLists.txt	Updates MODULE_ESTATE_elecstate_energy test linkage/sources for new symbols.

Comments suppressed due to low confidence (2)

source/source_estate/elecstate_energy.cpp:365

• This Makov–Payne block recomputes a full Hartree potential (H_Hartree_pw::v_hartree) and the correction every time cal_energies() is called. In the main SCF loop cal_energies(1) and cal_energies(2) are both called, so this work will be duplicated each iteration. Suggest computing/caching the correction once per potential update (or at least only in one of the two calls) and reusing it for the other energy functional.

    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<double> 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();
        }

source/source_estate/elecstate_energy.cpp:371

• The new Makov–Payne execution path (including acceptance of assume_isolated = mp/m-p/makov-payne and the actual correction applied to f_en.correction_el) is not covered by unit tests. Since source/source_estate/test/elecstate_energy_test.cpp already tests cal_energies(), please add at least one test that enables Makov–Payne and asserts that correction_el is non-zero/expected (and that the code path does not error for a supported cubic latname).

    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<double> 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;
    }

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[ErjieWu]

LGTM.
It seems it support only three kind of structures now and maybe we should add more in the future. For more precise values, could the code of the Ewald energy calculation module be reused? I'm not sure.
At the same time, this might be expanded into a separate class for electrostatic correction. More post-computation methods might be able to be implemented more uniformly through this class.

[mohanchen]

I accept current method, note that this feature needs further improvements.