delete omega in pelec, remove some deepks codes to setup_deepks.cpp

Author: mohanchen

source/source_esolver/esolver_double_xc.cpp

@@ -90,7 +90,6 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
 
     // 10) inititlize the charge density
     this->chr_base.allocate(PARAM.inp.nspin);
-    this->pelec_base->omega = ucell.omega;
 
     // 11) initialize the potential
     if (this->pelec_base->pot == nullptr)
@@ -116,7 +115,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 
     ESolver_KS_LCAO<TK,TR>::before_scf(ucell, istep);
 
-    this->pelec_base->omega = ucell.omega;
     //----------------------------------------------------------
     //! calculate D2 or D3 vdW
     //----------------------------------------------------------

source/source_esolver/esolver_fp.cpp

@@ -108,8 +108,6 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->locpp.init_vloc(ucell, this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
-        this->pelec->omega = ucell.omega;
-
         // perform symmetry analysis
         if (ModuleSymmetry::Symmetry::symm_flag == 1)
         {

source/source_esolver/esolver_ks_lcao.cpp

@@ -110,7 +110,6 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
 
     // 11) init charge density
     this->chr.allocate(inp.nspin);
-    this->pelec->omega = ucell.omega;
 
     // 12) init potentials
     if (this->pelec->pot == nullptr)
@@ -565,15 +564,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
-#ifdef __MLALGO
-    if (PARAM.inp.deepks_scf)
-    {
-        this->deepks.ld.dpks_cal_e_delta_band(dm_vec, this->kv.get_nks());
-        DeePKS_domain::update_dmr(this->kv.kvec_d, dm_vec, ucell, orb_, this->pv, this->gd, this->deepks.ld.dm_r);
-        estate->f_en.edeepks_scf = this->deepks.ld.E_delta - this->deepks.ld.e_delta_band;
-        estate->f_en.edeepks_delta = this->deepks.ld.E_delta;
-    }
-#endif
+    this->deepks.delta_e(this->kv, this->pv, this->gd, dm_vec, PARAM.inp);
 
     // 3) for delta spin
     if (PARAM.inp.sc_mag_switch)

source/source_esolver/esolver_of.cpp

@@ -273,7 +273,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
             {
                 // Here we initialize rho to be uniform,
                 // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements.
-                this->chr.rho[is][ibs] = this->nelec_[is] / this->pelec->omega;
+                this->chr.rho[is][ibs] = this->nelec_[is] / ucell.omega;
                 this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
             }
         }
@@ -587,7 +587,7 @@ void ESolver_OF::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
 {
     ModuleBase::matrix kinetic_stress_;
     kinetic_stress_.create(3, 3);
-    this->kedf_manager_->get_stress(this->pelec->omega,
+    this->kedf_manager_->get_stress(ucell.omega,
                                     this->chr.rho,
                                     this->pphi_,
                                     this->pw_rho,

source/source_esolver/esolver_of_tool.cpp

@@ -22,7 +22,6 @@ void ESolver_OF::init_elecstate(UnitCell& ucell)
         this->pelec = new elecstate::ElecState((Charge*)(&chr), this->pw_rho, pw_big);
         this->chr.allocate(PARAM.inp.nspin);
     }
-    this->pelec->omega = ucell.omega;
 
     delete this->pelec->pot;
     this->pelec->pot = new elecstate::Potential(this->pw_rhod,

source/source_estate/elecstate.h

@@ -114,7 +114,6 @@ class ElecState
     std::string classname = "elecstate";
 
     int iter = 0;                                  ///< scf iteration
-    double omega = 0.0;                            ///< volume
     Potential* pot = nullptr;                      ///< pointer to potential
     Charge* charge = nullptr;                      ///< pointer to charge density
     const K_Vectors* klist = nullptr;              ///< pointer to k points lists

source/source_estate/elecstate_energy.cpp

@@ -174,7 +174,7 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
     deband0 = deband_aux;
 #endif
 
-    deband0 *= this->omega / this->charge->rhopw->nxyz;
+    deband0 *= ucell.omega / this->charge->rhopw->nxyz;
 
     // \int rho(r) v_{exx}(r) dr = 2 E_{exx}[rho]
     deband0 -= 2 * this->f_en.exx; // Peize Lin add 2017-10-16
@@ -246,7 +246,7 @@ double ElecState::cal_delta_escf() const
 
     assert(this->charge->rhopw->nxyz > 0);
 
-    descf *= this->omega / this->charge->rhopw->nxyz;
+    descf *= this->charge->rhopw->omega / this->charge->rhopw->nxyz;
     return descf;
 }
 

source/source_estate/magnetism.cpp

@@ -23,8 +23,11 @@ void Magnetism::compute_mag(const double& omega,
 		const double* const * rho, 
 		double* nelec_spin)
 {
+    assert(omega>0.0);
     assert(nxyz>0);
 
+    const double fac = omega / nxyz;
+
     if (PARAM.inp.nspin==2)
     {
         this->tot_mag = 0.00;
@@ -40,8 +43,8 @@ void Magnetism::compute_mag(const double& omega,
         Parallel_Reduce::reduce_pool(this->tot_mag);
         Parallel_Reduce::reduce_pool(this->abs_mag);
 #endif
-        this->tot_mag *= omega / nxyz;
-        this->abs_mag *= omega / nxyz;
+        this->tot_mag *= fac;
+        this->abs_mag *= fac;
 
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Total magnetism (Bohr mag/cell)",this->tot_mag);
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Absolute magnetism (Bohr mag/cell)",this->abs_mag);
@@ -82,15 +85,15 @@ void Magnetism::compute_mag(const double& omega,
 #endif
 		for(int i=0;i<3;i++) 
 		{
-			this->tot_mag_nc[i] *= omega/ nxyz;
+			this->tot_mag_nc[i] *= fac;
             // mohan add 2025-06-21
 			if( std::abs(this->tot_mag_nc[i]) < 1.0e-16)
 			{
 				this->tot_mag_nc[i] = 0.0;
 			}
 		}
 
-		this->abs_mag *= omega/ nxyz;
+		this->abs_mag *= fac;
 
         // mohan update 2025-06-21
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Total magnetism (Bohr mag/cell)",

source/source_estate/setup_estate_pw.cpp

@@ -39,9 +39,6 @@ void elecstate::setup_estate_pw(UnitCell& ucell, // unitcell
         }
     }
 
-    //! Set the cell volume variable in pelec
-    pelec->omega = ucell.omega;
-
     //! Inititlize the charge density.
     chr.allocate(inp.nspin);
 

source/source_lcao/setup_deepks.cpp

@@ -61,6 +61,26 @@ void Setup_DeePKS<TK>::before_runner(const UnitCell& ucell, // unitcell
 #endif
 }
 
+template <typename TK>
+void Setup_DeePKS<TK>::delta_e(
+        const K_Vectors &kv,
+	    const Parallel_Orbitals &pv, // parallel orbitals
+		const Grid_Driver &gd,
+		const std::vector<std::vector<TK>>& dm_vec,
+		const Input_para &inp)
+{
+#ifdef __MLALGO
+    if (inp.deepks_scf)
+    {
+        this->ld.dpks_cal_e_delta_band(dm_vec, kv.get_nks());
+        DeePKS_domain::update_dmr(kv.kvec_d, dm_vec, ucell, orb_, pv, gd, this->ld.dm_r);
+        estate->f_en.edeepks_scf = this->ld.E_delta - this->ld.e_delta_band;
+        estate->f_en.edeepks_delta = this->ld.E_delta;
+    }
+#endif
+}
+
+
 
 template <typename TK>
 void Setup_DeePKS<TK>::write_forces(