Refactor: split pseudopot_cell_vl and pseudopot_cell_vnl

YuLiu98 · YuLiu98 · commit cddace7c489b · 2024-12-19T12:38:57.000+08:00
diff --git a/source/module_esolver/cell_dependency.h b/source/module_esolver/cell_dependency.h
diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp
@@ -285,7 +285,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
             this->pw_rhod->collect_uniqgg();
         }
 
-        this->p_locpp->init_vloc(ucell,this->pw_rhod);
+        this->locpp.init_vloc(ucell, this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
         this->pelec->omega = ucell.omega;
diff --git a/source/module_esolver/esolver_fp.h b/source/module_esolver/esolver_fp.h
@@ -7,6 +7,7 @@
 #include "module_elecstate/elecstate.h"
 #include "module_elecstate/module_charge/charge_extra.h"
 #include "module_hamilt_general/module_surchem/surchem.h"
+#include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 #include <fstream>
@@ -40,8 +41,12 @@ namespace ModuleESolver
         //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
         virtual void after_scf(UnitCell& ucell, const int istep);
 
-        //! Electronic states
-        elecstate::ElecState* pelec = nullptr;
+        //! ------------------------------------------------------------------------------
+        //! These pointers will be deleted in the free_pointers() function every ion step.
+        //! ------------------------------------------------------------------------------
+        elecstate::ElecState* pelec = nullptr; ///< Electronic states
+
+        //! ------------------------------------------------------------------------------
 
         //! Electorn charge density
         Charge chr;
@@ -54,8 +59,8 @@ namespace ModuleESolver
 
         ModulePW::PW_Basis* pw_rho;
 
-        //! pointer to pseudopotential
-        pseudopot_cell_vl* p_locpp = nullptr;
+        //! pointer to local pseudopotential
+        pseudopot_cell_vl locpp;
 
         /**
          * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense'
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
@@ -69,7 +69,6 @@ ESolver_KS<T, Device>::ESolver_KS()
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
     this->ppcell.cell_factor = PARAM.inp.cell_factor;
-    this->p_locpp = &this->ppcell;
 }
 
 //------------------------------------------------------------------------------
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -199,7 +199,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 
     // 8) initialize ppcell
-    this->ppcell.init_vloc(ucell, this->pw_rho);
+    this->locpp.init_vloc(ucell, this->pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     // 9) inititlize the charge density
@@ -212,7 +212,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         this->pelec->pot = new elecstate::Potential(this->pw_rhod,
                                                     this->pw_rho,
                                                     &ucell,
-                                                    &(this->ppcell.vloc),
+                                                    &(this->locpp.vloc),
                                                     &(this->sf),
                                                     &(this->solvent),
                                                     &(this->pelec->f_en.etxc),
@@ -317,7 +317,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        orb_,
                        force,
                        this->scs,
-                       this->ppcell,
+                       this->locpp,
                        this->sf,
                        this->kv,
                        this->pw_rho,
@@ -465,7 +465,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                     this->solvent,
                                     *this->pw_rho,
                                     *this->pw_rhod,
-                                    this->ppcell.vloc,
+                                    this->locpp.vloc,
                                     *this->pelec->charge,
                                     this->GG,
                                     this->GK,
@@ -493,7 +493,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                             this->solvent,
                                             *this->pw_rho,
                                             *this->pw_rhod,
-                                            this->ppcell.vloc,
+                                            this->locpp.vloc,
                                             *this->pelec->charge,
                                             this->GG,
                                             this->GK,
diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp
@@ -251,7 +251,7 @@ namespace ModuleESolver
                                 *this->pw_wfc,
                                 *this->pw_rho,
                                 *this->pw_rhod,
-                                this->ppcell.vloc,
+                                this->locpp.vloc,
                                 *this->pelec->charge,
                                 this->kv,
                                 this->pelec->wg
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
@@ -179,7 +179,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
         this->pelec->pot = new elecstate::Potential(this->pw_rhod,
                                                     this->pw_rho,
                                                     &ucell,
-                                                    &this->ppcell.vloc,
+                                                    &this->locpp.vloc,
                                                     &(this->sf),
                                                     &(this->solvent),
                                                     &(this->pelec->f_en.etxc),
@@ -206,14 +206,12 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
         delete this->psi;
     }
 
-    //! init pseudopotential
-    this->ppcell.init(ucell,&this->sf, this->pw_wfc);
-
     //! initalize local pseudopotential
-    this->ppcell.init_vloc(ucell,this->pw_rhod);
+    this->locpp.init_vloc(ucell, this->pw_rhod);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     //! Initalize non-local pseudopotential
+    this->ppcell.init(ucell, &this->sf, this->pw_wfc);
     this->ppcell.init_vnl(ucell, this->pw_rhod);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL");
 
@@ -308,7 +306,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     //! calculate the total local pseudopotential in real space
-    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc);
+    this->pelec->init_scf(istep, ucell, this->sf.strucFac, this->locpp.numeric, ucell.symm, (void*)this->pw_wfc);
 
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
@@ -780,7 +778,7 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                  &ucell.symm,
                  &this->sf,
                  this->solvent,
-                 &this->ppcell,
+                 &this->locpp,
                  &this->ppcell,
                  &this->kv,
                  this->pw_wfc,
@@ -802,6 +800,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(UnitCell& ucell, ModuleBase::matrix& s
                            : reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->kspw_psi);
     ss.cal_stress(stress,
                   ucell,
+                  this->locpp,
                   this->ppcell,
                   this->pw_rhod,
                   &ucell.symm,
diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp
@@ -23,7 +23,6 @@ ESolver_OF::ESolver_OF()
 {
     this->classname = "ESolver_OF";
     this->task_ = new char[60];
-    this->p_locpp = &this->locpp;
 }
 
 ESolver_OF::~ESolver_OF()
diff --git a/source/module_esolver/esolver_of.h b/source/module_esolver/esolver_of.h
@@ -30,8 +30,6 @@ class ESolver_OF : public ESolver_FP
     virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override;
 
     virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) override;
-  protected:
-    pseudopot_cell_vl locpp;
 
   private:
     // ======================= variables ==========================
diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp
@@ -233,6 +233,7 @@ void ESolver_SDFT_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix&
                     &this->sf,
                     &this->kv,
                     this->pw_wfc,
+                    this->locpp,
                     this->ppcell,
                     ucell,
                     *this->kspw_psi,
@@ -253,6 +254,7 @@ void ESolver_SDFT_PW<T, Device>::cal_stress(UnitCell& ucell, ModuleBase::matrix&
                   *this->kspw_psi,
                   this->stowf,
                   this->pelec->charge,
+                  &this->locpp,
                   &this->ppcell,
                   ucell);
 }
diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp
@@ -260,7 +260,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 #endif // __EXX
 
-    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm);
+    this->pelec->init_scf(istep, ucell, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
@@ -371,7 +371,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
         rdmft_solver.update_ion(ucell,
                                 *(this->pw_rho),
-                                this->ppcell.vloc,
+                                this->locpp.vloc,
                                 this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08
     }
 
diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp
@@ -249,7 +249,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     // pelec should be initialized before these calculations
-    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm);
+    this->pelec->init_scf(istep, ucell, this->sf.strucFac, this->locpp.numeric, ucell.symm);
     // self consistent calculations for electronic ground state
     if (cal_type == "get_pchg")
     {
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -46,7 +46,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                                           const LCAO_Orbitals& orb,
                                           ModuleBase::matrix& fcs,
                                           ModuleBase::matrix& scs,
-                                          const pseudopot_cell_vnl& nlpp,
+                                          const pseudopot_cell_vl& locpp,
                                           const Structure_Factor& sf,
                                           const K_Vectors& kv,
                                           ModulePW::PW_Basis* rhopw,
@@ -106,7 +106,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                              pelec->vnew_exist,
                              pelec->charge,
                              rhopw,
-                             nlpp,
+                             locpp,
                              sf);
     }
 
@@ -151,7 +151,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                               pelec->f_en.etxc,
                               pelec->charge,
                               rhopw,
-                              nlpp,
+                              locpp,
                               sf);
     }
 
@@ -280,7 +280,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
     if (PARAM.inp.imp_sol && isforce)
     {
         fsol.create(nat, 3);
-        solvent.cal_force_sol(ucell, rhopw, nlpp.vloc, fsol);
+        solvent.cal_force_sol(ucell, rhopw, locpp.vloc, fsol);
     }
 
     //! atomic forces from DFT+U (Quxin version)
@@ -840,15 +840,15 @@ void Force_Stress_LCAO<T>::calForcePwPart(const UnitCell& ucell,
                                           const bool vnew_exist,
                                           const Charge* const chr,
                                           ModulePW::PW_Basis* rhopw,
-                                          const pseudopot_cell_vnl& nlpp,
+                                          const pseudopot_cell_vl& locpp,
                                           const Structure_Factor& sf)
 {
     ModuleBase::TITLE("Force_Stress_LCAO", "calForcePwPart");
     //--------------------------------------------------------
     // local pseudopotential force:
     // use charge density; plane wave; local pseudopotential;
     //--------------------------------------------------------
-    f_pw.cal_force_loc(ucell, fvl_dvl, rhopw, nlpp.vloc, chr);
+    f_pw.cal_force_loc(ucell, fvl_dvl, rhopw, locpp.vloc, chr);
     //--------------------------------------------------------
     // ewald force: use plane wave only.
     //--------------------------------------------------------
@@ -857,11 +857,11 @@ void Force_Stress_LCAO<T>::calForcePwPart(const UnitCell& ucell,
     //--------------------------------------------------------
     // force due to core correlation.
     //--------------------------------------------------------
-    f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, GlobalC::ucell);
+    f_pw.cal_force_cc(fcc, rhopw, chr, locpp.numeric, GlobalC::ucell);
     //--------------------------------------------------------
     // force due to self-consistent charge.
     //--------------------------------------------------------
-    f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, nlpp.numeric, ucell);
+    f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, locpp.numeric, ucell);
     return;
 }
 
@@ -984,15 +984,15 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
                                            const double& etxc,
                                            const Charge* const chr,
                                            ModulePW::PW_Basis* rhopw,
-                                           const pseudopot_cell_vnl& nlpp,
+                                           const pseudopot_cell_vl& locpp,
                                            const Structure_Factor& sf)
 {
     ModuleBase::TITLE("Force_Stress_LCAO", "calStressPwPart");
     //--------------------------------------------------------
     // local pseudopotential stress:
     // use charge density; plane wave; local pseudopotential;
     //--------------------------------------------------------
-    sc_pw.stress_loc(ucell, sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr);
+    sc_pw.stress_loc(ucell, sigmadvl, rhopw, locpp.vloc, &sf, 0, chr);
 
     //--------------------------------------------------------
     // hartree term
@@ -1007,7 +1007,7 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
     //--------------------------------------------------------
     // stress due to core correlation.
     //--------------------------------------------------------
-    sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, nlpp.numeric, chr);
+    sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, locpp.numeric, chr);
 
     //--------------------------------------------------------
     // stress due to self-consistent charge.
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
@@ -44,7 +44,7 @@ class Force_Stress_LCAO
                         const LCAO_Orbitals& orb,
                         ModuleBase::matrix& fcs,
                         ModuleBase::matrix& scs,
-                        const pseudopot_cell_vnl& nlpp,
+                        const pseudopot_cell_vl& locpp,
                         const Structure_Factor& sf,
                         const K_Vectors& kv,
                         ModulePW::PW_Basis* rhopw,
@@ -76,7 +76,7 @@ class Force_Stress_LCAO
                         const bool vnew_exist,
                         const Charge* const chr,
                         ModulePW::PW_Basis* rhopw,
-                        const pseudopot_cell_vnl& nlpp,
+                        const pseudopot_cell_vl& locpp,
                         const Structure_Factor& sf);
 
     void integral_part(const bool isGammaOnly,
@@ -114,7 +114,7 @@ class Force_Stress_LCAO
                          const double& etxc,
                          const Charge* const chr,
                          ModulePW::PW_Basis* rhopw,
-                         const pseudopot_cell_vnl& nlpp,
+                         const pseudopot_cell_vl& locpp,
                          const Structure_Factor& sf);
 
     static double force_invalid_threshold_ev;
diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.h
@@ -2,6 +2,7 @@
 #define OF_STRESS_PW_H
 
 #include "module_elecstate/elecstate.h"
+#include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/stress_func.h"
 
 class OF_Stress_PW : public Stress_Func<double>
diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h
@@ -1,7 +1,6 @@
 #ifndef VNL_IN_PW_H
 #define VNL_IN_PW_H
 
-#include "VL_in_pw.h"
 #include "module_base/complexarray.h"
 #include "module_base/complexmatrix.h"
 #include "module_base/intarray.h"
@@ -18,7 +17,7 @@
 // Calculate the non-local pseudopotential in reciprocal
 // space using plane wave as basis set.
 //==========================================================
-class pseudopot_cell_vnl : public pseudopot_cell_vl
+class pseudopot_cell_vnl
 {
 
   public:
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h
@@ -9,6 +9,7 @@
 #include "module_cell/klist.h"
 #include "module_cell/module_symmetry/symmetry.h"
 #include "module_elecstate/elecstate.h"
+#include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/kernels/force_op.h"
 #include "module_hsolver/kernels/math_kernel_op.h"
 #include "module_psi/psi.h"
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.cpp
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h b/source/module_hamilt_pw/hamilt_pwdft/stress_pw.h
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_forces.cpp
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_forces.h b/source/module_hamilt_pw/hamilt_stodft/sto_forces.h
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.cpp
diff --git a/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h b/source/module_hamilt_pw/hamilt_stodft/sto_stress_pw.h

Original file line number	Diff line number	Diff line change
`@@ -285,7 +285,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)`
`285`	`285`	`this->pw_rhod->collect_uniqgg();`
`286`	`286`	`}`
`287`	`287`
`288`		`- this->p_locpp->init_vloc(ucell,this->pw_rhod);`
	`288`	`+ this->locpp.init_vloc(ucell, this->pw_rhod);`
`289`	`289`	`ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");`
`290`	`290`
`291`	`291`	`this->pelec->omega = ucell.omega;`
Original file line number	Diff line number	Diff line change
`@@ -69,7 +69,6 @@ ESolver_KS<T, Device>::ESolver_KS()`
`69`	`69`	`p_chgmix = new Charge_Mixing();`
`70`	`70`	`p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);`
`71`	`71`	`this->ppcell.cell_factor = PARAM.inp.cell_factor;`
`72`		`- this->p_locpp = &this->ppcell;`
`73`	`72`	`}`
`74`	`73`
`75`	`74`	`//------------------------------------------------------------------------------`
Original file line number	Diff line number	Diff line change
`@@ -23,7 +23,6 @@ ESolver_OF::ESolver_OF()`
`23`	`23`	`{`
`24`	`24`	`this->classname = "ESolver_OF";`
`25`	`25`	`this->task_ = new char[60];`
`26`		`- this->p_locpp = &this->locpp;`
`27`	`26`	`}`
`28`	`27`
`29`	`28`	`ESolver_OF::~ESolver_OF()`