Refactor: remove GlocalC::ucell in elecstate_energy_terms.cpp (#5630)

* Refactor: remove GlocalC::ucell in elecstate_energy_terms.cpp

* [pre-commit.ci lite] apply automatic fixes

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

Author: YuLiu98

source/Makefile.Objects

@@ -438,7 +438,6 @@ OBJS_SURCHEM=surchem.o\
     cal_totn.o\
     cal_vcav.o\
     cal_vel.o\
-    corrected_energy.o\
     minimize_cg.o\
     sol_force.o\
 

source/module_elecstate/elecstate_energy_terms.cpp

@@ -26,12 +26,12 @@ double ElecState::get_etot_gatefield()
 
 double ElecState::get_solvent_model_Ael()
 {
-    return GlobalC::solvent_model.cal_Ael(GlobalC::ucell, this->charge->nrxx, this->charge->nxyz);
+    return GlobalC::solvent_model.Ael;
 }
 
 double ElecState::get_solvent_model_Acav()
 {
-    return GlobalC::solvent_model.cal_Acav(GlobalC::ucell, this->charge->nxyz);
+    return GlobalC::solvent_model.Acav;
 }
 
 #ifdef __LCAO

source/module_hamilt_general/module_surchem/CMakeLists.txt

@@ -8,7 +8,6 @@ add_library(
     cal_totn.cpp
     cal_vcav.cpp
     cal_vel.cpp
-    corrected_energy.cpp
     minimize_cg.cpp
     sol_force.cpp
 )

source/module_hamilt_general/module_surchem/cal_vcav.cpp

@@ -10,18 +10,21 @@ void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW
     std::complex<double> *aux = new std::complex<double>[rho_basis->nmaxgr];
 
     // the formula is : rho(r)^prime = \int iG * rho(G)e^{iGr} dG
-    for (int ig = 0; ig < rho_basis->npw; ig++)
+    for (int ig = 0; ig < rho_basis->npw; ig++) {
         gdrtmpg[ig] = rhog[ig];
+}
     for(int i = 0 ; i < 3 ; ++i)
     {
         // calculate the charge density gradient in reciprocal space.
-        for (int ig = 0; ig < rho_basis->npw; ig++)
+        for (int ig = 0; ig < rho_basis->npw; ig++) {
             aux[ig] = gdrtmpg[ig] * pow(rho_basis->gcar[ig][i], 2);
+}
         // bring the gdr from G --> R
         rho_basis->recip2real(aux, aux);
         // remember to multily 2pi/a0, which belongs to G vectors.
-        for (int ir = 0; ir < rho_basis->nrxx; ir++)
+        for (int ir = 0; ir < rho_basis->nrxx; ir++) {
             lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
+}
     }
 
     delete[] gdrtmpg;
@@ -116,6 +119,8 @@ void surchem::createcavity(const UnitCell& ucell,
     //-------------------------------------------------------------
     // cavitation energy
     //-------------------------------------------------------------
+    this->Acav = PARAM.inp.tau * qs * ucell.omega / rho_basis->nxyz;
+    Parallel_Reduce::reduce_pool(this->Acav);
 
     // double Ael = cal_Acav(ucell, pwb);
 

source/module_hamilt_general/module_surchem/cal_vel.cpp

@@ -109,7 +109,13 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
     }
 
     // calculate Ael
-    // double Ael = cal_Ael(cell, pwb);
+    this->Ael = 0.0;
+    for (int ir = 0; ir < rho_basis->nrxx; ir++)
+    {
+        this->Ael -= TOTN_real[ir] * delta_phi[ir];
+    }
+    Parallel_Reduce::reduce_pool(this->Ael);
+    this->Ael *= cell.omega / rho_basis->nxyz;
 
     // the 2nd item of tmp_Vel
     eps_pot(PS_TOTN_real, Sol_phi, rho_basis, epsilon, epspot);

source/module_hamilt_general/module_surchem/corrected_energy.cpp

@@ -24,6 +24,9 @@ class surchem
     ModuleBase::matrix Vel;
     double qs;
 
+    double Acav;
+    double Ael;
+
     // get atom info
     atom_in GetAtom;