Electrostatic correction of Makov-Payne type

source/Makefile.Objects

@@ -234,6 +234,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\

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

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 <cmath>
@@ -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<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;
+    }
+
 #ifdef __MLALGO
     this->f_en.ml_exx = this->pot->get_ml_exx_energy();
 #endif

source/source_estate/fp_energy.cpp

@@ -16,23 +16,23 @@ 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;
 }
 
 /// @brief calculate etot_harris
 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;
 }
 
 /// @brief set all energies to zero
 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;

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