Update initialization of pw_wfc in esolver_ks.cpp, also some old gint has been deleted, E_bandgap has been changed to E_gap(k), struct efermi has been changed to Efermi (deepmodeling#6617)

* add update_pot in source_estate

* split update_pot in rt-TDDFT into two functions, one is the original update_pot, which is replaced by update_pot function in estate, the other is save2 function used in after_iter

* fix bugs

* update

* update esolver_ks

* small updates

* small bug fixed

* update E_bandgap to E_gap(k)

* delete rho_restart file

* change setup_parameters to print_parameters

* move setup_pw.cpp to setup_pwwfc.cpp in module_pwdft

* update esolver, small things

* delete some old_gint codes in esolver

* fix bugs, update outputs

* fix bug in tests

* change efermi to Efermi in fp_energy.h

* fix bugs

* fix exd and exc in ctrl_runner_lcao

* fix cmake

* fix bug

Author: mohanchen

source/Makefile.Objects

@@ -276,7 +276,6 @@ OBJS_ESOLVER=esolver.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
     pw_others.o\
-    pw_setup.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\
@@ -549,7 +548,7 @@ OBJS_IO=input_conv.o\
     cal_test.o\
     write_dos_pw.o\
     nscf_fermi_surf.o\
-    nscf_band.o\
+    write_bands.o\
     cal_dos.o\
     cal_pdos_gamma.o\
     cal_pdos_multik.o\
@@ -582,6 +581,7 @@ OBJS_IO=input_conv.o\
     write_elecstat_pot.o\
     write_elf.o\
     write_dipole.o\
+    write_init.o\
     td_current_io.o\
     write_libxc_r.o\
     output_log.o\
@@ -727,6 +727,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     fp_energy.o\
     setup_pot.o\
     setup_pwrho.o\
+    setup_pwwfc.o\
     forces.o\
     forces_us.o\
     forces_nl.o\

source/source_esolver/CMakeLists.txt

@@ -14,7 +14,6 @@ list(APPEND objects
     esolver_of_interface.cpp
     esolver_of_tool.cpp
     pw_others.cpp
-    pw_setup.cpp
 )
 if(ENABLE_LCAO)
   list(APPEND objects

source/source_esolver/esolver_fp.cpp

@@ -13,10 +13,11 @@
 #include "source_io/print_info.h"
 #include "source_io/rhog_io.h"
 #include "source_io/module_parameter/parameter.h"
-#include "source_io/ctrl_output_fp.h"
 
 #include "source_pw/module_pwdft/setup_pwrho.h" // mohan 20251005
 #include "source_hamilt/module_xc/xc_functional.h" // mohan 20251005
+#include "source_io/ctrl_output_fp.h"
+#include "source_io/write_init.h" // write_chg_init, write_pot_init
 
 namespace ModuleESolver
 {
@@ -145,64 +146,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     //! set direction of magnetism, used in non-collinear case 
     elecstate::cal_ux(ucell);
 
-    //! output the initial charge density
-    const int nspin = PARAM.inp.nspin;
-    if (PARAM.inp.out_chg[0] == 2)
-    {
-        for (int is = 0; is < nspin; is++)
-        {
-            std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "chg";
-
-			if(nspin==1)
-			{
-				ss << "ini.cube";
-			}
-			else if(nspin==2 || nspin==4)
-			{
-				ss << "s" << is + 1 << "ini.cube";
-			}
-
-            ModuleIO::write_vdata_palgrid(this->Pgrid,
-                                          this->chr.rho[is],
-                                          is,
-                                          nspin,
-                                          istep,
-                                          ss.str(),
-                                          this->pelec->eferm.ef,
-                                          &(ucell));
-        }
-    }
-
-    //! output total local potential of the initial charge density
-    if (PARAM.inp.out_pot == 3)
-    {
-        for (int is = 0; is < nspin; is++)
-        {
-            std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "pot";
-
-			if(nspin==1)
-			{
-				ss << "ini.cube";
-			}
-			else if(nspin==2 || nspin==4)
-			{
-				ss << "s" << is + 1 << "ini.cube";
-			}
-
-            ModuleIO::write_vdata_palgrid(this->Pgrid,
-                                          this->pelec->pot->get_effective_v(is),
-                                          is,
-                                          nspin,
-                                          istep,
-                                          ss.str(),
-                                          0.0, // efermi
-                                          &(ucell),
-                                          11, // precsion
-                                          0); // out_fermi
-        }
-    }
+    //! output the initial charge density and potential
+    ModuleIO::write_chg_init(ucell, this->Pgrid, this->chr, this->pelec->eferm, istep, PARAM.inp);
+//    ModuleIO::write_pot_init(ucell, this->Pgrid, this->pelec, istep, PARAM.inp); 
 
     return;
 }

source/source_esolver/esolver_gets.cpp

@@ -43,7 +43,7 @@ void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp)
     this->kv.set(ucell, ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
-    ModuleIO::setup_parameters(ucell, this->kv);
+    ModuleIO::print_parameters(ucell, this->kv, inp);
 
     // 2) init ElecState
     // autoset nbands in ElecState, it should before basis_init (for Psi 2d division)

source/source_esolver/esolver_ks.cpp

@@ -1,30 +1,19 @@
 #include "esolver_ks.h"
-#include "pw_setup.h" // setup plane wave
 
-#include "source_base/timer.h"
-#include "source_base/global_variable.h"
-#include "source_pw/module_pwdft/global.h"
-#include "source_io/module_parameter/parameter.h"
-#include "source_lcao/module_dftu/dftu.h"
-
-#include "source_cell/cal_atoms_info.h"
-#include "source_estate/elecstate_print.h"
-#include "source_hamilt/module_xc/xc_functional.h"
-#include "source_hsolver/hsolver.h"
-#include "source_io/cube_io.h"
-
-// for NSCF calculations of band structures
-#include "source_io/nscf_band.h"
-// for output log information
-#include "source_io/output_log.h"
-#include "source_io/print_info.h"
-#include "source_io/write_eig_occ.h"
 // for jason output information
 #include "source_io/json_output/init_info.h"
 #include "source_io/json_output/output_info.h"
 
 #include "source_estate/update_pot.h" // mohan add 20251016
 #include "source_estate/module_charge/chgmixing.h" // mohan add 20251018
+#include "source_pw/module_pwdft/setup_pwwfc.h" // mohan add 20251018
+#include "source_hsolver/hsolver.h"
+#include "source_io/write_eig_occ.h"
+#include "source_io/write_bands.h"
+#include "source_hamilt/module_xc/xc_functional.h"
+#include "source_io/output_log.h" // use write_head
+#include "source_estate/elecstate_print.h" // print_etot
+#include "source_io/print_info.h" // print_parameters
 
 namespace ModuleESolver
 {
@@ -40,10 +29,12 @@ ESolver_KS<T, Device>::~ESolver_KS()
 	// do not add any codes in this deconstructor funcion
 	//****************************************************
 	delete this->psi;
-    delete this->pw_wfc;
     delete this->p_hamilt;
     delete this->p_chgmix;
     this->ppcell.release_memory();
+    
+    // mohan add 2025-10-18, should be put int clean() function
+    pw::teardown_pwwfc(this->pw_wfc);
 }
 
 
@@ -59,74 +50,53 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     classname = "ESolver_KS";
     basisname = "";
 
-    scf_thr = inp.scf_thr;
-    scf_ene_thr = inp.scf_ene_thr;
-    maxniter = inp.scf_nmax;
-    niter = maxniter;
-    drho = 0.0;
-
-    std::string fft_device = inp.device;
+    this->scf_thr = inp.scf_thr;
+    this->scf_ene_thr = inp.scf_ene_thr;
+    this->maxniter = inp.scf_nmax;
+    this->niter = maxniter;
+    this->drho = 0.0;
 
-    //! 3) setup pw_wfc
-    // currently LCAO doesn't support GPU acceleration of FFT
-    if(inp.basis_type == "lcao")
-    {
-        fft_device = "cpu";
-    }
-    std::string fft_precision = inp.precision;
-#ifdef __ENABLE_FLOAT_FFTW
-    if (inp.cal_cond && inp.esolver_type == "sdft")
-    {
-        fft_precision = "mixing";
-    }
-#endif
-
-    pw_wfc = new ModulePW::PW_Basis_K_Big(fft_device, fft_precision);
-    ModulePW::PW_Basis_K_Big* tmp = static_cast<ModulePW::PW_Basis_K_Big*>(pw_wfc);
-
-    tmp->setbxyz(inp.bx, inp.by, inp.bz);
+    // cell_factor
+    this->ppcell.cell_factor = inp.cell_factor;
 
-    //! 4) setup charge mixing
+    //! 3) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
 
-    // cell_factor
-    this->ppcell.cell_factor = inp.cell_factor;
-
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    //! 5) setup Exc for the first element '0', because all elements have same exc 
+    //! 4) setup Exc for the first element '0', because all elements have same exc 
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
 
-    //! 6) setup the charge mixing parameters
+    //! 5) setup the charge mixing parameters
     p_chgmix->set_mixing(inp.mixing_mode, inp.mixing_beta, inp.mixing_ndim,
       inp.mixing_gg0, inp.mixing_tau, inp.mixing_beta_mag, inp.mixing_gg0_mag,
       inp.mixing_gg0_min, inp.mixing_angle, inp.mixing_dmr, ucell.omega, ucell.tpiba);
 
     p_chgmix->init_mixing();
 
-    //! 7) symmetry analysis should be performed every time the cell is changed
+    //! 6) symmetry analysis should be performed every time the cell is changed
     if (ModuleSymmetry::Symmetry::symm_flag == 1)
     {
         ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
 
-    //! 8) Setup the k points according to symmetry.
+    //! 7) Setup the k points according to symmetry.
     this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
-    //! 9) print information
-    ModuleIO::setup_parameters(ucell, this->kv);
+    //! 8) print information
+    ModuleIO::print_parameters(ucell, this->kv, inp);
 
-    //! 10) setup plane wave for electronic wave functions
-    ModuleESolver::pw_setup(inp, ucell, *this->pw_rho, this->kv, *this->pw_wfc);
+    //! 9) setup plane wave for electronic wave functions
+    pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);
 
-    //! 11) parallel of FFT grid 
+    //! 10) parallel of FFT grid 
 	Pgrid.init(this->pw_rhod->nx, this->pw_rhod->ny, this->pw_rhod->nz,
 			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
 
-    //! 12) calculate the structure factor
+    //! 11) calculate the structure factor
     this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod);
 }
 
@@ -303,13 +273,13 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
     }
 
     // 1.2) print out eigenvalues and occupations
-	if(iter % PARAM.inp.out_freq_elec == 0)
-	{
-		if (PARAM.inp.out_band[0] || iter == PARAM.inp.scf_nmax || conv_esolver)
+    if (PARAM.inp.out_band[0])
+    {
+		if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
 		{
 			ModuleIO::write_eig_iter(this->pelec->ekb,this->pelec->wg,*this->pelec->klist);
 		}
-	}
+    }
 
     // 2.1) compute magnetization, only for spin==2
     ucell.magnet.compute_mag(ucell.omega, this->chr.nrxx, this->chr.nxyz, this->chr.rho,
@@ -388,32 +358,9 @@ void ESolver_KS<T, Device>::after_scf(UnitCell& ucell, const int istep, const bo
     // 3) write eigenvalues and occupations to eig_occ.txt
     ModuleIO::write_eig_file(this->pelec->ekb, this->pelec->wg, this->kv, istep);
 
-    // 3) write band information to band.txt
-    if (PARAM.inp.out_band[0])
-    {
-        const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1;
-        for (int is = 0; is < nspin0; is++)
-        {
-            std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "band";
-
-            if(nspin0==1)
-            {
-                // do nothing
-            }
-            else if(nspin0==2)
-            {
-                ss << "s" << is + 1;
-            }
-
-            ss << ".txt";
+    // 4) write band information to band.txt
+    ModuleIO::write_bands(PARAM.inp, this->pelec->ekb, this->kv);
 
-            const double eshift = 0.0;
-            ModuleIO::nscf_band(is, ss.str(), PARAM.inp.nbands,
-                                eshift, PARAM.inp.out_band[1], // precision
-                                this->pelec->ekb, this->kv);
-        }
-    }
 }
 
 template <typename T, typename Device>

source/source_esolver/esolver_ks_lcao.cpp

@@ -294,12 +294,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
 		  this->gd, this->psi, this->chr, hamilt_lcao,
           this->two_center_bundle_, this->GG, this->GK,
           this->orb_, this->pw_rho, this->pw_rhod,
-          this->sf, this->locpp.vloc, 
-#ifdef __EXX
-				this->exx_nao.exd,
-				this->exx_nao.exc,
-#endif
-				this->solvent);
+          this->sf, this->locpp.vloc, this->exx_nao, this->solvent);
 
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners");
 }
@@ -313,9 +308,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     ESolver_KS<TK>::iter_init(ucell, istep, iter);
 
     // cast pointers
-
 	auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-
 	if(!estate)
 	{
 		ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_init","pelec does not exist");

source/source_esolver/esolver_ks_lcao.h

@@ -77,10 +77,8 @@ class ESolver_KS_LCAO : public ESolver_KS<TK>
     //! Grid integration: used to store some basic information
     Grid_Technique GridT;
 
-#ifndef __OLD_GINT
     //! GintInfo: used to store some basic infomation about module_gint
     std::unique_ptr<ModuleGint::GintInfo> gint_info_;
-#endif
 
     //! NAO orbitals: two-center integrations
     TwoCenterBundle two_center_bundle_;