From 32b673dceb79c4aceb70f31fdeb8ea451e847176 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Tue, 6 Jan 2026 15:41:12 -0700 Subject: [PATCH 1/3] sol refactor --- .../source_estate/module_pot/pot_surchem.hpp | 19 ++- .../module_surchem/H_correction_pw.cpp | 32 ++-- .../module_surchem/cal_pseudo.cpp | 4 +- .../source_hamilt/module_surchem/cal_vcav.cpp | 40 +++-- .../source_hamilt/module_surchem/cal_vel.cpp | 19 ++- .../module_surchem/minimize_cg.cpp | 145 +++++++++--------- .../source_hamilt/module_surchem/surchem.cpp | 5 +- source/source_hamilt/module_surchem/surchem.h | 48 +++--- .../module_surchem/test/cal_pseudo_test.cpp | 5 +- .../module_surchem/test/cal_vcav_test.cpp | 14 +- .../module_surchem/test/cal_vel_test.cpp | 11 +- 11 files changed, 193 insertions(+), 149 deletions(-) diff --git a/source/source_estate/module_pot/pot_surchem.hpp b/source/source_estate/module_pot/pot_surchem.hpp index 3ccb4960dc..dc06e60fd4 100644 --- a/source/source_estate/module_pot/pot_surchem.hpp +++ b/source/source_estate/module_pot/pot_surchem.hpp @@ -31,6 +31,7 @@ class PotSurChem : public PotBase } } + // Passing an explicit output matrix makes the lifetime and allocation explicit and avoids hidden allocations. void cal_v_eff(const Charge*const chg, const UnitCell*const ucell, ModuleBase::matrix& v_eff) override { if (!this->allocated) @@ -38,14 +39,16 @@ class PotSurChem : public PotBase this->surchem_->allocate(this->rho_basis_->nrxx, v_eff.nr); this->allocated = true; } - - v_eff += this->surchem_->v_correction(*ucell, - *chg->pgrid, - const_cast(this->rho_basis_), - v_eff.nr, - chg->rho, - this->vlocal, - this->structure_factors_); + ModuleBase::matrix v_sol_correction(v_eff.nr, this->rho_basis_->nrxx); + this->surchem_->v_correction(*ucell, + *chg->pgrid, + const_cast(this->rho_basis_), + v_eff.nr, + chg->rho, + this->vlocal, + this->structure_factors_, + v_sol_correction); + v_eff += v_sol_correction; } private: diff --git a/source/source_hamilt/module_surchem/H_correction_pw.cpp b/source/source_hamilt/module_surchem/H_correction_pw.cpp index 06bb2eea25..289f00d07e 100644 --- a/source/source_hamilt/module_surchem/H_correction_pw.cpp +++ b/source/source_hamilt/module_surchem/H_correction_pw.cpp @@ -5,14 +5,15 @@ #include "source_base/timer.h" #include "source_hamilt/module_xc/xc_functional.h" #include "surchem.h" - -ModuleBase::matrix surchem::v_correction(const UnitCell& cell, - const Parallel_Grid& pgrid, - const ModulePW::PW_Basis* rho_basis, - const int& nspin, - const double* const* const rho, - const double* vlocal, - Structure_Factor* sf) +// Changing the interface to use an explicit output parameter clarifies lifetime management and avoids hidden allocations. +void surchem::v_correction(const UnitCell& cell, + const Parallel_Grid& pgrid, + const ModulePW::PW_Basis* rho_basis, + const int& nspin, + const double* const* const rho, + const double* vlocal, + Structure_Factor* sf, + ModuleBase::matrix& v) { ModuleBase::TITLE("surchem", "v_cor"); ModuleBase::timer::tick("surchem", "v_cor"); @@ -46,10 +47,15 @@ ModuleBase::matrix surchem::v_correction(const UnitCell& cell, cal_pseudo(cell, pgrid, rho_basis, porter_g, ps_totn, sf); - ModuleBase::matrix v(nspin, rho_basis->nrxx); + // ModuleBase::matrix v(nspin, rho_basis->nrxx); + if (v.nr != nspin || v.nc != rho_basis->nrxx) + { + v.create(nspin, rho_basis->nrxx); + } + ModuleBase::GlobalFunc::ZEROS(v.c, nspin * rho_basis->nrxx); - v += cal_vel(cell, rho_basis, total_n, ps_totn, nspin); - v += cal_vcav(cell, rho_basis, ps_totn, nspin); + cal_vel(cell, rho_basis, total_n, ps_totn, nspin, v); + cal_vcav(cell, rho_basis, ps_totn, nspin, v); delete[] porter; delete[] porter_g; @@ -58,5 +64,5 @@ ModuleBase::matrix surchem::v_correction(const UnitCell& cell, delete[] total_n; ModuleBase::timer::tick("surchem", "v_cor"); - return v; -} + return; +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/cal_pseudo.cpp b/source/source_hamilt/module_surchem/cal_pseudo.cpp index 9f707654d0..9c3a8be429 100644 --- a/source/source_hamilt/module_surchem/cal_pseudo.cpp +++ b/source/source_hamilt/module_surchem/cal_pseudo.cpp @@ -9,8 +9,8 @@ void surchem::gauss_charge(const UnitCell& cell, std::complex* N, Structure_Factor* sf) { - sf->setup(&cell, pgrid, rho_basis); // this is strange, should be removed to other places, mohan add 2025-11-04 - + //sf->setup(&cell, pgrid, rho_basis); // this is strange, should be removed to other places, mohan add 2025-11-04 + //sf here was only needed in the test. const int ig0 = rho_basis->ig_gge0; // G=0 index for (int it = 0; it < cell.ntype; it++) { diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index 617da1376e..ecf213cfcd 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -16,25 +16,28 @@ void lapl_rho(const double& tpiba2, // the formula is : rho(r)^prime = \int iG * rho(G)e^{iGr} dG 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. + ModuleBase::GlobalFunc::ZEROS(aux, rho_basis->nmaxgr); + 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++) { lapn[ir] -= aux[ir].real() * tpiba2; -} - + } } delete[] gdrtmpg; delete[] aux; - return; } // calculates first derivative of the shape function in realspace @@ -66,6 +69,7 @@ void surchem::createcavity(const UnitCell& ucell, { ModuleBase::Vector3 *nablan = new ModuleBase::Vector3[rho_basis->nrxx]; ModuleBase::GlobalFunc::ZEROS(nablan, rho_basis->nrxx); + double *nablan_2 = new double[rho_basis->nrxx]; double *sqrt_nablan_2 = new double[rho_basis->nrxx]; double *lapn = new double[rho_basis->nrxx]; @@ -77,14 +81,14 @@ void surchem::createcavity(const UnitCell& ucell, // nabla n XC_Functional::grad_rho(ps_totn, nablan, rho_basis, ucell.tpiba); - // |\nabla n |^2 = nablan_2 + // |\nabla n |^2 = nablan_2 for (int ir = 0; ir < rho_basis->nrxx; ir++) { nablan_2[ir] = pow(nablan[ir].x, 2) + pow(nablan[ir].y, 2) + pow(nablan[ir].z, 2); } // Laplacian of n - lapl_rho(ucell.tpiba2,ps_totn, lapn, rho_basis); + lapl_rho(ucell.tpiba2, ps_totn, lapn, rho_basis); //------------------------------------------------------------- // add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2) @@ -162,10 +166,13 @@ void surchem::createcavity(const UnitCell& ucell, delete[] ggn; } -ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, - const ModulePW::PW_Basis* rho_basis, - std::complex* ps_totn, - int nspin) +//The interface is changed to use an explicit output parameter to +//clarify lifetime management and avoid hidden allocations. +void surchem::cal_vcav(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis, + std::complex* ps_totn, + int nspin, + ModuleBase::matrix& v) { ModuleBase::TITLE("surchem", "cal_vcav"); ModuleBase::timer::tick("surchem", "cal_vcav"); @@ -175,12 +182,12 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, createcavity(ucell, rho_basis, ps_totn, tmp_Vcav); - ModuleBase::GlobalFunc::ZEROS(Vcav.c, nspin * rho_basis->nrxx); if (nspin == 4) { for (int ir = 0; ir < rho_basis->nrxx; ir++) { - Vcav(0, ir) += tmp_Vcav[ir]; + Vcav(0, ir) = tmp_Vcav[ir]; + v(0, ir) += Vcav(0, ir); } } else @@ -189,12 +196,13 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, { for (int ir = 0; ir < rho_basis->nrxx; ir++) { - Vcav(is, ir) += tmp_Vcav[ir]; + Vcav(is, ir) = tmp_Vcav[ir]; + v(is, ir) += Vcav(is, ir); } } } delete[] tmp_Vcav; ModuleBase::timer::tick("surchem", "cal_vcav"); - return Vcav; -} + return; +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index 7790dd3341..83b72966a2 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -54,11 +54,14 @@ void eps_pot(const double* PS_TOTN_real, delete[] phisq; } -ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, - const ModulePW::PW_Basis* rho_basis, - std::complex* TOTN, - std::complex* PS_TOTN, - int nspin) +//The interface is changed to use an explicit output parameter to +//clarify lifetime management and avoid hidden allocations. +void surchem::cal_vel(const UnitCell& cell, + const ModulePW::PW_Basis* rho_basis, + std::complex* TOTN, + std::complex* PS_TOTN, + int nspin, + ModuleBase::matrix& v) { ModuleBase::TITLE("surchem", "cal_vel"); ModuleBase::timer::tick("surchem", "cal_vel"); @@ -134,6 +137,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, for (int ir = 0; ir < rho_basis->nrxx; ir++) { Vel(0, ir) += tmp_Vel[ir]; + v(0, ir) += Vel(0, ir); } } else @@ -143,6 +147,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, for (int ir = 0; ir < rho_basis->nrxx; ir++) { Vel(is, ir) += tmp_Vel[ir]; + v(is, ir) += Vel(is, ir); } } } @@ -158,5 +163,5 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, delete[] phi_tilda_R0; ModuleBase::timer::tick("surchem", "cal_vel"); - return Vel; -} + return; +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/minimize_cg.cpp b/source/source_hamilt/module_surchem/minimize_cg.cpp index fb9fe63f73..186c949b79 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -15,9 +15,8 @@ void surchem::minimize_cg(const UnitCell& ucell, double rinvLr = 0; // r * r double r2 = 0; - // precond loop parameter - int i = 0; ModuleBase::GlobalFunc::ZEROS(phi, rho_basis->npw); + // malloc vectors in G space std::complex *resid = new std::complex[rho_basis->npw]; std::complex *z = new std::complex[rho_basis->npw]; @@ -29,7 +28,17 @@ void surchem::minimize_cg(const UnitCell& ucell, std::complex *gradphi_y = new std::complex[rho_basis->npw]; std::complex *gradphi_z = new std::complex[rho_basis->npw]; - std::complex *phi_work = new std::complex[rho_basis->npw]; + // Removed unused phi_work allocation + // std::complex *phi_work = new std::complex[rho_basis->npw]; + + // ========================================================== + // PRE-ALLOCATION FOR LEPS2 (Avoids allocation inside loop) + // ========================================================== + ModuleBase::Vector3 *aux_grad_phi = new ModuleBase::Vector3[rho_basis->nrxx]; + std::complex *aux_grad_grad_phi_G = new std::complex[rho_basis->npw]; + ModuleBase::Vector3 *aux_tmp_vector3 = new ModuleBase::Vector3[rho_basis->nrxx]; + double *aux_lp_real = new double[rho_basis->nrxx]; + double *aux_grad_grad_phi_real = new double[rho_basis->nrxx]; ModuleBase::GlobalFunc::ZEROS(resid, rho_basis->npw); ModuleBase::GlobalFunc::ZEROS(z, rho_basis->npw); @@ -41,8 +50,6 @@ void surchem::minimize_cg(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(gradphi_y, rho_basis->npw); ModuleBase::GlobalFunc::ZEROS(gradphi_z, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(phi_work, rho_basis->npw); - int count = 0; double gg = 0; @@ -65,7 +72,9 @@ void surchem::minimize_cg(const UnitCell& ucell, } // call leps to calculate div ( epsilon * grad ) phi - Leps2(ucell, rho_basis, phi, d_eps, gradphi_x, gradphi_y, gradphi_z, phi_work, lp); + // Updated Leps2 call with new buffers + Leps2(ucell, rho_basis, phi, d_eps, gradphi_x, gradphi_y, gradphi_z, lp, + aux_grad_phi, aux_grad_grad_phi_G, aux_tmp_vector3, aux_lp_real, aux_grad_grad_phi_real); // the residue // r = A*phi + (chtot + N) @@ -85,8 +94,6 @@ void surchem::minimize_cg(const UnitCell& ucell, rinvLr = ModuleBase::GlobalFunc::ddot_real(rho_basis->npw, resid, z); r2 = ModuleBase::GlobalFunc::ddot_real(rho_basis->npw, resid, resid); - double r20 = r2; - // copy for (int ig = 0; ig < rho_basis->npw; ig++) { @@ -103,9 +110,10 @@ void surchem::minimize_cg(const UnitCell& ucell, break; } - Leps2(ucell, rho_basis, d, d_eps, gradphi_x, gradphi_y, gradphi_z, phi_work, lp); + // Updated Leps2 call inside loop + Leps2(ucell, rho_basis, d, d_eps, gradphi_x, gradphi_y, gradphi_z, lp, + aux_grad_phi, aux_grad_grad_phi_G, aux_tmp_vector3, aux_lp_real, aux_grad_grad_phi_real); - // cout <<"lp after leps"<npw, d, lp); // update phi @@ -149,7 +157,7 @@ void surchem::minimize_cg(const UnitCell& ucell, // output: num of cg loop ncgsol = count; - // comment test res + // CLEANUP delete[] resid; delete[] z; delete[] lp; @@ -158,87 +166,82 @@ void surchem::minimize_cg(const UnitCell& ucell, delete[] gradphi_x; delete[] gradphi_y; delete[] gradphi_z; - delete[] phi_work; + // delete[] phi_work; // Removed + + // Clean up auxiliary buffers + delete[] aux_grad_phi; + delete[] aux_grad_grad_phi_G; + delete[] aux_tmp_vector3; + delete[] aux_lp_real; + delete[] aux_grad_grad_phi_real; } +// avoid creating large temporary matrices inside its iteration loop void surchem::Leps2(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, std::complex* phi, - double* epsilon, // epsilon from shapefunc, dim=nrxx - std::complex* gradphi_x, // dim=ngmc + double* epsilon, // epsilon from shapefunc, dim=nrxx + std::complex* gradphi_x, std::complex* gradphi_y, std::complex* gradphi_z, - std::complex* phi_work, - std::complex* lp) + std::complex* lp, + ModuleBase::Vector3* grad_phi_R, // size: nrxx + std::complex* aux_G, // size: npw + ModuleBase::Vector3* tmp_vector3, // size: nrxx) + double* lp_real, // size: nrxx + double* aux_R) // size: nrxx { - ModuleBase::Vector3 *grad_phi = new ModuleBase::Vector3[rho_basis->nrxx]; - XC_Functional::grad_rho(phi, grad_phi, rho_basis, ucell.tpiba); - - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_phi[ir].x *= epsilon[ir]; - grad_phi[ir].y *= epsilon[ir]; - grad_phi[ir].z *= epsilon[ir]; - } - std::vector lp_real(rho_basis->nrxx,0); - ModuleBase::GlobalFunc::ZEROS(lp, rho_basis->npw); + XC_Functional::grad_rho(phi, grad_phi_R, rho_basis, ucell.tpiba); - std::vector grad_grad_phi(rho_basis->nrxx,0); - std::complex *grad_grad_phi_G = new std::complex[rho_basis->npw]; - ModuleBase::Vector3 *tmp_vector3 = new ModuleBase::Vector3[rho_basis->nrxx]; - // x - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_G, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_grad_phi[ir] = grad_phi[ir].x; - } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].x; + grad_phi_R[ir].x *= epsilon[ir]; + grad_phi_R[ir].y *= epsilon[ir]; + grad_phi_R[ir].z *= epsilon[ir]; } - // y - grad_grad_phi.assign(grad_grad_phi.size(),0.0); - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_G, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_grad_phi[ir] = grad_phi[ir].y; + + ModuleBase::GlobalFunc::ZEROS(lp_real, rho_basis->nrxx); + + // 1. R -> G + for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].x; + rho_basis->real2recip(aux_R, gradphi_x); // + + + for(int ig=0; ignpw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_x[ig] * rho_basis->gcar[ig][0]; // 0 = x } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].y; + rho_basis->recip2real(aux_G, aux_R); + for(int ir=0; irnrxx; ir++) { + lp_real[ir] += aux_R[ir] * ucell.tpiba; } - // z - grad_grad_phi.assign(grad_grad_phi.size(),0.0); - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_G, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_grad_phi[ir] = grad_phi[ir].z; + + for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].y; + rho_basis->real2recip(aux_R, gradphi_y); + + for(int ig=0; ignpw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_y[ig] * rho_basis->gcar[ig][1]; // 1 = y } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].z; + rho_basis->recip2real(aux_G, aux_R); + for(int ir=0; irnrxx; ir++) { + lp_real[ir] += aux_R[ir] * ucell.tpiba; } - rho_basis->real2recip(lp_real.data(), lp); + for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].z; + rho_basis->real2recip(aux_R, gradphi_z); + + for(int ig=0; ignpw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_z[ig] * rho_basis->gcar[ig][2]; // 2 = z + } + rho_basis->recip2real(aux_G, aux_R); + for(int ir=0; irnrxx; ir++) { + lp_real[ir] += aux_R[ir] * ucell.tpiba; + } - delete[] grad_phi; - std::vector().swap(lp_real); - std::vector().swap(grad_grad_phi); - delete[] grad_grad_phi_G; - delete[] tmp_vector3; -} + rho_basis->real2recip(lp_real, lp); +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/surchem.cpp b/source/source_hamilt/module_surchem/surchem.cpp index 60db451b41..1f2a66ef86 100644 --- a/source/source_hamilt/module_surchem/surchem.cpp +++ b/source/source_hamilt/module_surchem/surchem.cpp @@ -50,9 +50,12 @@ void surchem::clear() this->TOTN_real = nullptr; this->delta_phi = nullptr; this->epspot = nullptr; + + this->Vcav.create(0, 0); + this->Vel.create(0, 0); } surchem::~surchem() { this->clear(); -} +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/surchem.h b/source/source_hamilt/module_surchem/surchem.h index e2358449c1..383c72b7d0 100644 --- a/source/source_hamilt/module_surchem/surchem.h +++ b/source/source_hamilt/module_surchem/surchem.h @@ -62,16 +62,18 @@ class surchem const std::complex* PS_TOTN, double* vwork); - ModuleBase::matrix cal_vcav(const UnitCell& ucell, - const ModulePW::PW_Basis* rho_basis, - std::complex* PS_TOTN, - int nspin); + void cal_vcav(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis, + std::complex* PS_TOTN, + int nspin, + ModuleBase::matrix& v); - ModuleBase::matrix cal_vel(const UnitCell& cell, - const ModulePW::PW_Basis* rho_basis, - std::complex* TOTN, - std::complex* PS_TOTN, - int nspin); + void cal_vel(const UnitCell& cell, + const ModulePW::PW_Basis* rho_basis, + std::complex* TOTN, + std::complex* PS_TOTN, + int nspin, + ModuleBase::matrix& v); double cal_Ael(const UnitCell& cell, const int& nrxx, // num. of real space grids on current core @@ -99,16 +101,22 @@ class surchem std::complex* gradphi_x, // dim=ngmc std::complex* gradphi_y, std::complex* gradphi_z, - std::complex* phi_work, - std::complex* lp); - - ModuleBase::matrix v_correction(const UnitCell& cell, - const Parallel_Grid& pgrid, - const ModulePW::PW_Basis* rho_basis, - const int& nspin, - const double* const* const rho, - const double* vlocal, - Structure_Factor* sf); + std::complex* lp, + // New buffers + ModuleBase::Vector3* grad_phi, + std::complex* grad_grad_phi_G, + ModuleBase::Vector3* tmp_vector3, + double* lp_real, + double* grad_grad_phi_real); + + void v_correction(const UnitCell& cell, + const Parallel_Grid& pgrid, + const ModulePW::PW_Basis* rho_basis, + const int& nspin, + const double* const* const rho, + const double* vlocal, + Structure_Factor* sf, + ModuleBase::matrix& v); void test_V_to_N(ModuleBase::matrix& v, const UnitCell& cell, @@ -134,4 +142,4 @@ class surchem private: }; -#endif +#endif \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/test/cal_pseudo_test.cpp b/source/source_hamilt/module_surchem/test/cal_pseudo_test.cpp index e60b192719..2bd8025f7d 100644 --- a/source/source_hamilt/module_surchem/test/cal_pseudo_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_pseudo_test.cpp @@ -76,6 +76,8 @@ TEST_F(cal_pseudo_test, gauss_charge) Structure_Factor sf; sf.nbspline = -1; + sf.setup(&ucell, pgrid, &pwtest); + solvent_model.gauss_charge(ucell, pgrid, &pwtest, N, &sf); EXPECT_NEAR(N[14].real(), 0.002, 1e-9); @@ -126,6 +128,7 @@ TEST_F(cal_pseudo_test, cal_pseudo) Structure_Factor sf; sf.nbspline = -1; + sf.setup(&ucell, pgrid, &pwtest); // sf.setup is moved to here std::complex* Porter_g = new std::complex[npw]; ModuleBase::GlobalFunc::ZEROS(Porter_g, npw); @@ -160,4 +163,4 @@ int main(int argc, char** argv) #endif return result; -} +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp b/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp index 65fe839aa2..9a0af0432d 100644 --- a/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp @@ -253,12 +253,14 @@ TEST_F(cal_vcav_test, cal_vcav) int nspin = 2; solvent_model.Vcav.create(nspin, nrxx); - solvent_model.cal_vcav(ucell, &pwtest, PS_TOTN, nspin); + ModuleBase::matrix v_res(nspin, nrxx); + ModuleBase::GlobalFunc::ZEROS(v_res.c, nspin * nrxx); + solvent_model.cal_vcav(ucell, &pwtest, PS_TOTN, nspin, v_res); - EXPECT_NEAR(solvent_model.Vcav(0, 0), 4.8556305312, 1e-10); - EXPECT_NEAR(solvent_model.Vcav(0, 1), -2.1006480538, 1e-10); - EXPECT_NEAR(solvent_model.Vcav(1, 0), 4.8556305312, 1e-10); - EXPECT_NEAR(solvent_model.Vcav(1, 1), -2.1006480538, 1e-10); + EXPECT_NEAR(v_res(0, 0), 4.8556305312, 1e-10); + EXPECT_NEAR(v_res(0, 1), -2.1006480538, 1e-10); + EXPECT_NEAR(v_res(1, 0), 4.8556305312, 1e-10); + EXPECT_NEAR(v_res(1, 1), -2.1006480538, 1e-10); delete[] PS_TOTN; } @@ -279,4 +281,4 @@ int main(int argc, char** argv) #endif return result; -} +} \ No newline at end of file diff --git a/source/source_hamilt/module_surchem/test/cal_vel_test.cpp b/source/source_hamilt/module_surchem/test/cal_vel_test.cpp index 60d925b3fb..818e52dbe6 100644 --- a/source/source_hamilt/module_surchem/test/cal_vel_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_vel_test.cpp @@ -220,10 +220,13 @@ TEST_F(cal_vel_test, cal_vel) solvent_model.TOTN_real = new double[nrxx]; solvent_model.delta_phi = new double[nrxx]; - solvent_model.cal_vel(ucell, &pwtest, TOTN, PS_TOTN, nspin); + ModuleBase::matrix v_res(nspin, nrxx); + ModuleBase::GlobalFunc::ZEROS(v_res.c, nspin * nrxx); - EXPECT_NEAR(solvent_model.Vel(0, 0), 0.0532168705, 1e-10); - EXPECT_NEAR(solvent_model.Vel(0, 1), 0.0447818244, 1e-10); + solvent_model.cal_vel(ucell, &pwtest, TOTN, PS_TOTN, nspin, v_res); + + EXPECT_NEAR(v_res(0, 0), 0.0532168705, 1e-10); + EXPECT_NEAR(v_res(0, 1), 0.0447818244, 1e-10); delete[] PS_TOTN; delete[] TOTN; @@ -245,4 +248,4 @@ int main(int argc, char** argv) #endif return result; -} +} \ No newline at end of file From 9256a90f889522d092a5291f6ee144f8395bb4b6 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Tue, 6 Jan 2026 17:00:40 -0700 Subject: [PATCH 2/3] remove unused variable --- .../module_surchem/minimize_cg.cpp | 41 +++++++------------ source/source_hamilt/module_surchem/surchem.h | 7 +--- 2 files changed, 16 insertions(+), 32 deletions(-) diff --git a/source/source_hamilt/module_surchem/minimize_cg.cpp b/source/source_hamilt/module_surchem/minimize_cg.cpp index 186c949b79..451a8fc233 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -24,9 +24,7 @@ void surchem::minimize_cg(const UnitCell& ucell, std::complex *gsqu = new std::complex[rho_basis->npw]; std::complex *d = new std::complex[rho_basis->npw]; - std::complex *gradphi_x = new std::complex[rho_basis->npw]; - std::complex *gradphi_y = new std::complex[rho_basis->npw]; - std::complex *gradphi_z = new std::complex[rho_basis->npw]; + std::complex *gradphi_G_work = new std::complex[rho_basis->npw]; // Removed unused phi_work allocation // std::complex *phi_work = new std::complex[rho_basis->npw]; @@ -36,7 +34,6 @@ void surchem::minimize_cg(const UnitCell& ucell, // ========================================================== ModuleBase::Vector3 *aux_grad_phi = new ModuleBase::Vector3[rho_basis->nrxx]; std::complex *aux_grad_grad_phi_G = new std::complex[rho_basis->npw]; - ModuleBase::Vector3 *aux_tmp_vector3 = new ModuleBase::Vector3[rho_basis->nrxx]; double *aux_lp_real = new double[rho_basis->nrxx]; double *aux_grad_grad_phi_real = new double[rho_basis->nrxx]; @@ -45,10 +42,7 @@ void surchem::minimize_cg(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(lp, rho_basis->npw); ModuleBase::GlobalFunc::ZEROS(gsqu, rho_basis->npw); ModuleBase::GlobalFunc::ZEROS(d, rho_basis->npw); - - ModuleBase::GlobalFunc::ZEROS(gradphi_x, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(gradphi_y, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(gradphi_z, rho_basis->npw); + ModuleBase::GlobalFunc::ZEROS(gradphi_G_work, rho_basis->npw); int count = 0; double gg = 0; @@ -73,8 +67,8 @@ void surchem::minimize_cg(const UnitCell& ucell, // call leps to calculate div ( epsilon * grad ) phi // Updated Leps2 call with new buffers - Leps2(ucell, rho_basis, phi, d_eps, gradphi_x, gradphi_y, gradphi_z, lp, - aux_grad_phi, aux_grad_grad_phi_G, aux_tmp_vector3, aux_lp_real, aux_grad_grad_phi_real); + Leps2(ucell, rho_basis, phi, d_eps, gradphi_G_work, lp, + aux_grad_phi, aux_grad_grad_phi_G, aux_lp_real, aux_grad_grad_phi_real); // the residue // r = A*phi + (chtot + N) @@ -111,8 +105,8 @@ void surchem::minimize_cg(const UnitCell& ucell, } // Updated Leps2 call inside loop - Leps2(ucell, rho_basis, d, d_eps, gradphi_x, gradphi_y, gradphi_z, lp, - aux_grad_phi, aux_grad_grad_phi_G, aux_tmp_vector3, aux_lp_real, aux_grad_grad_phi_real); + Leps2(ucell, rho_basis, d, d_eps, gradphi_G_work, lp, + aux_grad_phi, aux_grad_grad_phi_G, aux_lp_real, aux_grad_grad_phi_real); // calculate alpha alpha = -rinvLr / ModuleBase::GlobalFunc::ddot_real(rho_basis->npw, d, lp); @@ -163,15 +157,11 @@ void surchem::minimize_cg(const UnitCell& ucell, delete[] lp; delete[] gsqu; delete[] d; - delete[] gradphi_x; - delete[] gradphi_y; - delete[] gradphi_z; - // delete[] phi_work; // Removed + delete[] gradphi_G_work; // Clean up auxiliary buffers delete[] aux_grad_phi; delete[] aux_grad_grad_phi_G; - delete[] aux_tmp_vector3; delete[] aux_lp_real; delete[] aux_grad_grad_phi_real; } @@ -181,13 +171,10 @@ void surchem::Leps2(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, std::complex* phi, double* epsilon, // epsilon from shapefunc, dim=nrxx - std::complex* gradphi_x, - std::complex* gradphi_y, - std::complex* gradphi_z, + std::complex* gradphi_G_work, std::complex* lp, ModuleBase::Vector3* grad_phi_R, // size: nrxx std::complex* aux_G, // size: npw - ModuleBase::Vector3* tmp_vector3, // size: nrxx) double* lp_real, // size: nrxx double* aux_R) // size: nrxx { @@ -207,11 +194,11 @@ void surchem::Leps2(const UnitCell& ucell, // 1. R -> G for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].x; - rho_basis->real2recip(aux_R, gradphi_x); // + rho_basis->real2recip(aux_R, gradphi_G_work); // for(int ig=0; ignpw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_x[ig] * rho_basis->gcar[ig][0]; // 0 = x + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_G_work[ig] * rho_basis->gcar[ig][0]; // 0 = x } rho_basis->recip2real(aux_G, aux_R); for(int ir=0; irnrxx; ir++) { @@ -220,10 +207,10 @@ void surchem::Leps2(const UnitCell& ucell, for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].y; - rho_basis->real2recip(aux_R, gradphi_y); + rho_basis->real2recip(aux_R, gradphi_G_work); for(int ig=0; ignpw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_y[ig] * rho_basis->gcar[ig][1]; // 1 = y + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_G_work[ig] * rho_basis->gcar[ig][1]; // 1 = y } rho_basis->recip2real(aux_G, aux_R); for(int ir=0; irnrxx; ir++) { @@ -232,10 +219,10 @@ void surchem::Leps2(const UnitCell& ucell, for (int ir = 0; ir < rho_basis->nrxx; ir++) aux_R[ir] = grad_phi_R[ir].z; - rho_basis->real2recip(aux_R, gradphi_z); + rho_basis->real2recip(aux_R, gradphi_G_work); for(int ig=0; ignpw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_z[ig] * rho_basis->gcar[ig][2]; // 2 = z + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_G_work[ig] * rho_basis->gcar[ig][2]; // 2 = z } rho_basis->recip2real(aux_G, aux_R); for(int ir=0; irnrxx; ir++) { diff --git a/source/source_hamilt/module_surchem/surchem.h b/source/source_hamilt/module_surchem/surchem.h index 383c72b7d0..9c057aaffb 100644 --- a/source/source_hamilt/module_surchem/surchem.h +++ b/source/source_hamilt/module_surchem/surchem.h @@ -98,14 +98,11 @@ class surchem const ModulePW::PW_Basis* rho_basis, std::complex* phi, double* epsilon, // epsilon from shapefunc, dim=nrxx - std::complex* gradphi_x, // dim=ngmc - std::complex* gradphi_y, - std::complex* gradphi_z, + std::complex* gradphi_G_work, std::complex* lp, // New buffers - ModuleBase::Vector3* grad_phi, + ModuleBase::Vector3* grad_phi_R, std::complex* grad_grad_phi_G, - ModuleBase::Vector3* tmp_vector3, double* lp_real, double* grad_grad_phi_real); From c640326e1c209625e527bcd6732680964919dd30 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Tue, 6 Jan 2026 17:40:37 -0700 Subject: [PATCH 3/3] Clean up comments in cal_vcav.cpp --- source/source_hamilt/module_surchem/cal_vcav.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index ecf213cfcd..32f4db6f57 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -20,9 +20,7 @@ void lapl_rho(const double& tpiba2, for(int i = 0 ; i < 3 ; ++i) { - // calculate the charge density gradient in reciprocal space. - ModuleBase::GlobalFunc::ZEROS(aux, rho_basis->nmaxgr); - + // calculate the charge density gradient in reciprocal space. for (int ig = 0; ig < rho_basis->npw; ig++) { aux[ig] = gdrtmpg[ig] * pow(rho_basis->gcar[ig][i], 2); } @@ -205,4 +203,4 @@ void surchem::cal_vcav(const UnitCell& ucell, delete[] tmp_Vcav; ModuleBase::timer::tick("surchem", "cal_vcav"); return; -} \ No newline at end of file +}