From f474d9f6ab47c01943d023b5aefb9284dfcbe23c Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Sun, 4 Jan 2026 16:09:13 -0700 Subject: [PATCH 1/9] sol_refine --- .../source_estate/module_pot/pot_surchem.hpp | 18 +-- .../module_surchem/H_correction_pw.cpp | 28 ++-- .../module_surchem/cal_pseudo.cpp | 2 +- .../source_hamilt/module_surchem/cal_vcav.cpp | 13 +- .../source_hamilt/module_surchem/cal_vel.cpp | 7 +- .../module_surchem/minimize_cg.cpp | 122 +++++++++++------- source/source_hamilt/module_surchem/surchem.h | 46 ++++--- 7 files changed, 142 insertions(+), 94 deletions(-) diff --git a/source/source_estate/module_pot/pot_surchem.hpp b/source/source_estate/module_pot/pot_surchem.hpp index 3ccb4960dc..dd1c87d594 100644 --- a/source/source_estate/module_pot/pot_surchem.hpp +++ b/source/source_estate/module_pot/pot_surchem.hpp @@ -38,14 +38,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..aaabc23b14 100644 --- a/source/source_hamilt/module_surchem/H_correction_pw.cpp +++ b/source/source_hamilt/module_surchem/H_correction_pw.cpp @@ -6,13 +6,14 @@ #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) +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; } diff --git a/source/source_hamilt/module_surchem/cal_pseudo.cpp b/source/source_hamilt/module_surchem/cal_pseudo.cpp index 9f707654d0..b5855979d8 100644 --- a/source/source_hamilt/module_surchem/cal_pseudo.cpp +++ b/source/source_hamilt/module_surchem/cal_pseudo.cpp @@ -9,7 +9,7 @@ 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 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..2accbba9d0 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -162,10 +162,11 @@ 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) +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"); @@ -181,6 +182,7 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, for (int ir = 0; ir < rho_basis->nrxx; ir++) { Vcav(0, ir) += tmp_Vcav[ir]; + v(0, ir) += Vcav(0, ir); } } else @@ -190,11 +192,12 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell, for (int ir = 0; ir < rho_basis->nrxx; 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; } diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index 7790dd3341..c68491ffe4 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -58,7 +58,8 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, std::complex* TOTN, std::complex* PS_TOTN, - int nspin) + int nspin, + ModuleBase::matrix& v) { ModuleBase::TITLE("surchem", "cal_vel"); ModuleBase::timer::tick("surchem", "cal_vel"); @@ -134,6 +135,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 +145,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 +161,5 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell, delete[] phi_tilda_R0; ModuleBase::timer::tick("surchem", "cal_vel"); - return Vel; + return; } diff --git a/source/source_hamilt/module_surchem/minimize_cg.cpp b/source/source_hamilt/module_surchem/minimize_cg.cpp index fb9fe63f73..778fa14078 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -16,8 +16,10 @@ void surchem::minimize_cg(const UnitCell& ucell, // r * r double r2 = 0; // precond loop parameter - int i = 0; + // int i = 0; // Unused + 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 +31,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 +53,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 +75,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 +97,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 +113,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 +160,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 +169,102 @@ 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; } 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, + std::complex* gradphi_x, std::complex* gradphi_y, std::complex* gradphi_z, - std::complex* phi_work, - std::complex* lp) + std::complex* lp, + // New arguments for memory buffers + ModuleBase::Vector3* grad_phi, + std::complex* grad_grad_phi_G, + ModuleBase::Vector3* tmp_vector3, + double* lp_real, + double* grad_grad_phi_real) { - ModuleBase::Vector3 *grad_phi = new ModuleBase::Vector3[rho_basis->nrxx]; + // RESET BUFFERS at the start of call + // Previously these were "new" allocations, now we must clear them manually + ModuleBase::GlobalFunc::ZEROS(grad_phi, rho_basis->nrxx); + ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_G, rho_basis->npw); + ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); + ModuleBase::GlobalFunc::ZEROS(lp_real, rho_basis->nrxx); + // grad_grad_phi_real is overwritten by assignment, so ZEROS is optional but safe + ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); + ModuleBase::GlobalFunc::ZEROS(lp, rho_basis->npw); + // Calculate grad_rho XC_Functional::grad_rho(phi, grad_phi, rho_basis, ucell.tpiba); + // Multiply by epsilon 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); - 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); + // --- X Component --- for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_grad_phi[ir] = grad_phi[ir].x; + grad_grad_phi_real[ir] = grad_phi[ir].x; } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); + rho_basis->real2recip(grad_grad_phi_real, grad_grad_phi_G); + + // reuse tmp_vector3 + ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); 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; } - // 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); + // --- Y Component --- + ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); // Clear buffer for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_grad_phi[ir] = grad_phi[ir].y; + grad_grad_phi_real[ir] = grad_phi[ir].y; } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); + rho_basis->real2recip(grad_grad_phi_real, grad_grad_phi_G); + + ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); 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; } - // 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); + // --- Z Component --- + ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); // Clear buffer for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_grad_phi[ir] = grad_phi[ir].z; + grad_grad_phi_real[ir] = grad_phi[ir].z; } - rho_basis->real2recip(grad_grad_phi.data(), grad_grad_phi_G); + rho_basis->real2recip(grad_grad_phi_real, grad_grad_phi_G); + + ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); 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; } + // Final transfer to Reciprocal space + rho_basis->real2recip(lp_real, lp); - rho_basis->real2recip(lp_real.data(), lp); - - delete[] grad_phi; - std::vector().swap(lp_real); - std::vector().swap(grad_grad_phi); - - delete[] grad_grad_phi_G; - delete[] tmp_vector3; -} + // No delete[] calls here anymore! +} \ 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..7afcae58ab 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, From 474c40373bedc2a2fcfa4da20e9315430fc96e35 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Sun, 4 Jan 2026 16:19:36 -0700 Subject: [PATCH 2/9] fix void --- source/source_hamilt/module_surchem/cal_vel.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index c68491ffe4..eed83a5c8a 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -54,12 +54,12 @@ 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, - ModuleBase::matrix& v) +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"); From 4d60da8ee7b107c1220fb650e96e175b77d46372 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 01:18:20 -0700 Subject: [PATCH 3/9] mem_leak_try --- .../source_hamilt/module_surchem/cal_vcav.cpp | 132 ++++++++++------ .../module_surchem/minimize_cg.cpp | 142 ++++++++++-------- 2 files changed, 159 insertions(+), 115 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index 2accbba9d0..e5bfb90519 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -3,6 +3,36 @@ #include "source_io/module_parameter/parameter.h" #include "surchem.h" + +static void local_grad_rho(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis, + const std::complex* rho_G, // + ModuleBase::Vector3* grad_R, // + std::complex* aux_G, // + double* aux_R) // +{ + + for (int i = 0; i < 3; ++i) + { + // 1. iG * rho(G) + for (int ig = 0; ig < rho_basis->npw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; + } + + // 2. FFT: G -> R + + rho_basis->recip2real(aux_G, aux_R); + + // 3. 2pi/a + for (int ir = 0; ir < rho_basis->nrxx; ir++) { + grad_R[ir][i] = aux_R[ir] * ucell.tpiba; + } + } +} + +// --------------------------------------------------------- +// Helper 2: Laplacian +// --------------------------------------------------------- void lapl_rho(const double& tpiba2, const std::complex* rhog, double* lapn, @@ -11,37 +41,41 @@ void lapl_rho(const double& tpiba2, std::complex *gdrtmpg = new std::complex[rho_basis->npw]; ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx); - std::complex *aux = new std::complex[rho_basis->nmaxgr]; - // the formula is : rho(r)^prime = \int iG * rho(G)e^{iGr} dG + std::complex *aux = new std::complex[rho_basis->nrxx]; + + 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->nrxx); + 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 -// exp(-(log(n/n_c))^2 /(2 sigma^2)) /(sigma * sqrt(2*pi) )/n +// --------------------------------------------------------- +// Helper 3: Shape function +// --------------------------------------------------------- void shape_gradn(const std::complex* ps_totn, const ModulePW::PW_Basis* rho_basis, double* eprime) { - double *ps_totn_real = new double[rho_basis->nrxx]; ModuleBase::GlobalFunc::ZEROS(ps_totn_real, rho_basis->nrxx); rho_basis->recip2real(ps_totn, ps_totn_real); @@ -59,6 +93,9 @@ void shape_gradn(const std::complex* ps_totn, const ModulePW::PW_Basis* delete[] ps_totn_real; } +// --------------------------------------------------------- +// Create Cavity +// --------------------------------------------------------- void surchem::createcavity(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, const std::complex* ps_totn, @@ -66,6 +103,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]; @@ -74,22 +112,21 @@ void surchem::createcavity(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(sqrt_nablan_2, rho_basis->nrxx); ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx); - // nabla n - XC_Functional::grad_rho(ps_totn, nablan, rho_basis, ucell.tpiba); - // |\nabla n |^2 = nablan_2 + std::complex *aux_G = new std::complex[rho_basis->npw]; + double *aux_R = new double[rho_basis->nrxx]; + + local_grad_rho(ucell, rho_basis, ps_totn, nablan, aux_G, aux_R); + + // | \nabla n |^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); - //------------------------------------------------------------- - // add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2) - // to sqrt_nablan_2 - //------------------------------------------------------------- + lapl_rho(ucell.tpiba2, ps_totn, lapn, rho_basis); + double tmp = 0; double min = 1e-10; @@ -100,48 +137,40 @@ void surchem::createcavity(const UnitCell& ucell, sqrt_nablan_2[ir] = tmp; } - //------------------------------------------------------------- - // term1 = gamma*A / n, where - // gamma * A = exp(-(log(n/n_c))^2 /(2 sigma^2)) /(sigma * sqrt(2*pi) ) - //------------------------------------------------------------- + // -------------------------------------------------------- + // term1 = gamma*A / n + // -------------------------------------------------------- double *term1 = new double[rho_basis->nrxx]; shape_gradn(ps_totn, rho_basis, term1); - //------------------------------------------------------------- - // quantum surface area, integral of (gamma*A / n) * |\nabla n| - //=term1 * sqrt_nablan_2 - //------------------------------------------------------------- + // -------------------------------------------------------- + // quantum surface area calculation + // -------------------------------------------------------- qs = 0; - for (int ir = 0; ir < rho_basis->nrxx; ir++) { qs = qs + (term1[ir]) * (sqrt_nablan_2[ir]); - - // 1 / |nabla n| + // 1 / |nabla n| sqrt_nablan_2[ir] = 1 / std::max(sqrt_nablan_2[ir], min); } - //------------------------------------------------------------- - // cavitation energy - //------------------------------------------------------------- + // 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); - // packs the real array into a complex one - // to G space std::complex *inv_gn = new std::complex[rho_basis->npw]; rho_basis->real2recip(sqrt_nablan_2, inv_gn); - // \nabla(1 / |\nabla n|), ggn in real space + ModuleBase::Vector3 *ggn = new ModuleBase::Vector3[rho_basis->nrxx]; - XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba); + + + local_grad_rho(ucell, rho_basis, inv_gn, ggn, aux_G, aux_R); - //------------------------------------------------------------- - // add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav in real space - // and multiply by term1 = gamma*A/n in real space - //------------------------------------------------------------- + // -------------------------------------------------------- + // add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav + // -------------------------------------------------------- for (int ir = 0; ir < rho_basis->nrxx; ir++) { tmp = nablan[ir].x * ggn[ir].x + nablan[ir].y * ggn[ir].y + nablan[ir].z * ggn[ir].z; @@ -153,6 +182,9 @@ void surchem::createcavity(const UnitCell& ucell, vwork[ir] = vwork[ir] * term1[ir] * PARAM.inp.tau; } + + delete[] aux_G; + delete[] aux_R; delete[] nablan; delete[] nablan_2; delete[] sqrt_nablan_2; @@ -162,6 +194,9 @@ void surchem::createcavity(const UnitCell& ucell, delete[] ggn; } +// --------------------------------------------------------- +// Main Function: cal_vcav +// --------------------------------------------------------- void surchem::cal_vcav(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, std::complex* ps_totn, @@ -176,12 +211,11 @@ void 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); } } @@ -191,7 +225,7 @@ void 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); } } @@ -200,4 +234,4 @@ void surchem::cal_vcav(const UnitCell& ucell, delete[] tmp_Vcav; ModuleBase::timer::tick("surchem", "cal_vcav"); 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 778fa14078..b14ac263d7 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -16,7 +16,7 @@ void surchem::minimize_cg(const UnitCell& ucell, // r * r double r2 = 0; // precond loop parameter - // int i = 0; // Unused + ModuleBase::GlobalFunc::ZEROS(phi, rho_basis->npw); @@ -179,92 +179,102 @@ void surchem::minimize_cg(const UnitCell& ucell, delete[] aux_grad_grad_phi_real; } +void helper_grad_rho(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis, + const std::complex* rho_G, // G + ModuleBase::Vector3* grad_rho_R, // R + std::complex* buffer_G, + double* buffer_R) +{ + // 1.x, y, z + for (int i = 0; i < 3; ++i) + { + // 2. + + for (int ig = 0; ig < rho_basis->npw; ig++) + { + buffer_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; + } + + // 3. FFT: G -> R + + rho_basis->recip2real(buffer_G, buffer_R); + + + for (int ir = 0; ir < rho_basis->nrxx; ir++) + { + grad_rho_R[ir][i] = buffer_R[ir] * ucell.tpiba; + } + } +} + void surchem::Leps2(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, std::complex* phi, - double* epsilon, + double* epsilon, std::complex* gradphi_x, std::complex* gradphi_y, std::complex* gradphi_z, std::complex* lp, - // New arguments for memory buffers - ModuleBase::Vector3* grad_phi, - std::complex* grad_grad_phi_G, - ModuleBase::Vector3* tmp_vector3, - double* lp_real, - double* grad_grad_phi_real) + 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 { - // RESET BUFFERS at the start of call - // Previously these were "new" allocations, now we must clear them manually - ModuleBase::GlobalFunc::ZEROS(grad_phi, rho_basis->nrxx); - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_G, rho_basis->npw); - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - ModuleBase::GlobalFunc::ZEROS(lp_real, rho_basis->nrxx); - // grad_grad_phi_real is overwritten by assignment, so ZEROS is optional but safe - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); - ModuleBase::GlobalFunc::ZEROS(lp, rho_basis->npw); - // Calculate grad_rho - XC_Functional::grad_rho(phi, grad_phi, rho_basis, ucell.tpiba); - - // Multiply by epsilon - 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]; - } + helper_grad_rho(ucell, rho_basis, phi, grad_phi_R, aux_G, aux_R); + - // --- X Component --- for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_grad_phi_real[ir] = grad_phi[ir].x; + grad_phi_R[ir].x *= epsilon[ir]; + grad_phi_R[ir].y *= epsilon[ir]; + grad_phi_R[ir].z *= epsilon[ir]; } - rho_basis->real2recip(grad_grad_phi_real, grad_grad_phi_G); - - // reuse tmp_vector3 - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); + + + 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 ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].x; + + for(int ig=0; ignpw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_x[ig] * rho_basis->gcar[ig][0]; // 0 = x } - - // --- Y Component --- - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); // Clear buffer - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_grad_phi_real[ir] = grad_phi[ir].y; + rho_basis->recip2real(aux_G, aux_R); + for(int ir=0; irnrxx; ir++) { + lp_real[ir] += aux_R[ir] * ucell.tpiba; } - rho_basis->real2recip(grad_grad_phi_real, grad_grad_phi_G); - - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); + + + 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 ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].y; + for(int ig=0; ignpw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * gradphi_y[ig] * rho_basis->gcar[ig][1]; // 1 = y } - - // --- Z Component --- - ModuleBase::GlobalFunc::ZEROS(grad_grad_phi_real, rho_basis->nrxx); // Clear buffer - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_grad_phi_real[ir] = grad_phi[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(grad_grad_phi_real, grad_grad_phi_G); - - ModuleBase::GlobalFunc::ZEROS(tmp_vector3, rho_basis->nrxx); - XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis, ucell.tpiba); + + + 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 ir = 0; ir < rho_basis->nrxx; ir++) - { - lp_real[ir] += tmp_vector3[ir].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; } - // Final transfer to Reciprocal space - rho_basis->real2recip(lp_real, lp); - // No delete[] calls here anymore! + rho_basis->real2recip(lp_real, lp); } \ No newline at end of file From 896e63dddb1ac346bb0a537206161960f014270f Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 02:32:51 -0700 Subject: [PATCH 4/9] mem_leak --- source/source_estate/module_pot/efield.cpp | 2 +- .../source_hamilt/module_surchem/cal_vel.cpp | 65 +++++++++++++++++-- .../source_hamilt/module_surchem/surchem.cpp | 3 + 3 files changed, 65 insertions(+), 5 deletions(-) diff --git a/source/source_estate/module_pot/efield.cpp b/source/source_estate/module_pot/efield.cpp index c51b6262ac..34d24f7d49 100644 --- a/source/source_estate/module_pot/efield.cpp +++ b/source/source_estate/module_pot/efield.cpp @@ -227,7 +227,7 @@ double Efield::cal_induced_dipole(const UnitCell& cell, Parallel_Reduce::reduce_pool(induced_dipole); induced_dipole *= cell.lat0 / bmod * ModuleBase::FOUR_PI / rho_basis->nxyz; - + delete[] induced_rho; return induced_dipole; } diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index eed83a5c8a..cbb3c5e240 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -3,9 +3,33 @@ #include "source_io/module_parameter/parameter.h" #include "surchem.h" -void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis, double* eprime) + +static void local_grad_rho_vel(const UnitCell& ucell, + const ModulePW::PW_Basis* rho_basis, + const std::complex* rho_G, // G + ModuleBase::Vector3* grad_R, // R + std::complex* aux_G, // + double* aux_R) // { + for (int i = 0; i < 3; ++i) + { + // 1. iG * rho(G) + for (int ig = 0; ig < rho_basis->npw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; + } + // 2. FFT: G -> R + rho_basis->recip2real(aux_G, aux_R); + + // 3. 2pi/a + for (int ir = 0; ir < rho_basis->nrxx; ir++) { + grad_R[ir][i] = aux_R[ir] * ucell.tpiba; + } + } +} + +void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis, double* eprime) +{ double epr_c = 1.0 / sqrt(ModuleBase::TWO_PI) / PARAM.inp.sigma_k; double epr_z = 0; double min = 1e-10; @@ -16,6 +40,25 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis } } +static void local_grad_rho_vel_tpiba(const double tpiba, + const ModulePW::PW_Basis* rho_basis, + const std::complex* rho_G, + ModuleBase::Vector3* grad_R, + std::complex* aux_G, + double* aux_R) +{ + for (int i = 0; i < 3; ++i) + { + for (int ig = 0; ig < rho_basis->npw; ig++) { + aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; + } + rho_basis->recip2real(aux_G, aux_R); + for (int ir = 0; ir < rho_basis->nrxx; ir++) { + grad_R[ir][i] = aux_R[ir] * tpiba; + } + } +} + void eps_pot(const double* PS_TOTN_real, const double& tpiba, const std::complex* phi, @@ -36,8 +79,21 @@ void eps_pot(const double* PS_TOTN_real, ModuleBase::Vector3 *nabla_phi = new ModuleBase::Vector3[rho_basis->nrxx]; double *phisq = new double[rho_basis->nrxx]; - // nabla phi - XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba); + + + std::complex *aux_G = new std::complex[rho_basis->npw]; + double *aux_R = new double[rho_basis->nrxx]; + + + local_grad_rho_vel_tpiba(tpiba, rho_basis, phi, nabla_phi, aux_G, aux_R); + + + delete[] aux_G; + delete[] aux_R; + + // ======================================================== + // 修复结束 + // ======================================================== for (int ir = 0; ir < rho_basis->nrxx; ir++) { @@ -120,6 +176,7 @@ void surchem::cal_vel(const UnitCell& cell, this->Ael *= cell.omega / rho_basis->nxyz; // the 2nd item of tmp_Vel + // eps_pot internally calls shape_gradn and now the SAFE local_grad_rho_vel_tpiba eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot); for (int i = 0; i < rho_basis->nrxx; i++) @@ -162,4 +219,4 @@ void surchem::cal_vel(const UnitCell& cell, ModuleBase::timer::tick("surchem", "cal_vel"); return; -} +} \ 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..0b8067eb65 100644 --- a/source/source_hamilt/module_surchem/surchem.cpp +++ b/source/source_hamilt/module_surchem/surchem.cpp @@ -50,6 +50,9 @@ 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() From 3bd74b0c57d4db1199a701e1faa18be10ca3f2fe Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 02:41:37 -0700 Subject: [PATCH 5/9] clear --- source/source_hamilt/module_surchem/cal_vcav.cpp | 12 +++++------- source/source_hamilt/module_surchem/cal_vel.cpp | 12 +++++------- source/source_hamilt/module_surchem/minimize_cg.cpp | 10 ++-------- 3 files changed, 12 insertions(+), 22 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index e5bfb90519..4212a96816 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -6,10 +6,10 @@ static void local_grad_rho(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, // - ModuleBase::Vector3* grad_R, // - std::complex* aux_G, // - double* aux_R) // + const std::complex* rho_G, + ModuleBase::Vector3* grad_R, + std::complex* aux_G, + double* aux_R) { for (int i = 0; i < 3; ++i) @@ -30,9 +30,7 @@ static void local_grad_rho(const UnitCell& ucell, } } -// --------------------------------------------------------- -// Helper 2: Laplacian -// --------------------------------------------------------- + void lapl_rho(const double& tpiba2, const std::complex* rhog, double* lapn, diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index cbb3c5e240..5022296b3b 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -6,10 +6,10 @@ static void local_grad_rho_vel(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, // G - ModuleBase::Vector3* grad_R, // R - std::complex* aux_G, // - double* aux_R) // + const std::complex* rho_G, + ModuleBase::Vector3* grad_R, + std::complex* aux_G, + double* aux_R) { for (int i = 0; i < 3; ++i) { @@ -91,9 +91,7 @@ void eps_pot(const double* PS_TOTN_real, delete[] aux_G; delete[] aux_R; - // ======================================================== - // 修复结束 - // ======================================================== + for (int ir = 0; ir < rho_basis->nrxx; ir++) { diff --git a/source/source_hamilt/module_surchem/minimize_cg.cpp b/source/source_hamilt/module_surchem/minimize_cg.cpp index b14ac263d7..3581f861e3 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -8,14 +8,12 @@ void surchem::minimize_cg(const UnitCell& ucell, std::complex* phi, int& ncgsol) { - // parameters of CG method + double alpha = 0; double beta = 0; - // r * r' double rinvLr = 0; - // r * r double r2 = 0; - // precond loop parameter + ModuleBase::GlobalFunc::ZEROS(phi, rho_basis->npw); @@ -186,18 +184,14 @@ void helper_grad_rho(const UnitCell& ucell, std::complex* buffer_G, double* buffer_R) { - // 1.x, y, z for (int i = 0; i < 3; ++i) { - // 2. for (int ig = 0; ig < rho_basis->npw; ig++) { buffer_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; } - // 3. FFT: G -> R - rho_basis->recip2real(buffer_G, buffer_R); From dbd1a2b046b9fca2cfd8dcedace4bbc61ca64006 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 03:54:53 -0700 Subject: [PATCH 6/9] update test --- .../module_surchem/test/cal_vcav_test.cpp | 14 ++++++++------ .../module_surchem/test/cal_vel_test.cpp | 9 ++++++--- 2 files changed, 14 insertions(+), 9 deletions(-) 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..9163601080 100644 --- a/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp @@ -251,14 +251,16 @@ TEST_F(cal_vcav_test, cal_vcav) } int nspin = 2; - solvent_model.Vcav.create(nspin, nrxx); + // 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; } 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..1fe63e7fc9 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; From bbec623226504fece543085ed6cf6118c5b83e86 Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 05:10:26 -0700 Subject: [PATCH 7/9] fix test --- source/source_hamilt/module_surchem/test/cal_pseudo_test.cpp | 4 ++++ source/source_hamilt/module_surchem/test/cal_vcav_test.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) 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..bb0ca8b311 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); @@ -127,6 +129,8 @@ TEST_F(cal_pseudo_test, cal_pseudo) Structure_Factor sf; sf.nbspline = -1; + sf.setup(&ucell, pgrid, &pwtest); + std::complex* Porter_g = new std::complex[npw]; ModuleBase::GlobalFunc::ZEROS(Porter_g, npw); for (int i = 0; i < npw; i++) 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 9163601080..a864f0d152 100644 --- a/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp +++ b/source/source_hamilt/module_surchem/test/cal_vcav_test.cpp @@ -251,7 +251,7 @@ TEST_F(cal_vcav_test, cal_vcav) } int nspin = 2; - // solvent_model.Vcav.create(nspin, nrxx); + solvent_model.Vcav.create(nspin, nrxx); ModuleBase::matrix v_res(nspin, nrxx); ModuleBase::GlobalFunc::ZEROS(v_res.c, nspin * nrxx); From 4725616feb9ecd9e20e9f1c46ef30cd8939c053e Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 15:42:57 -0700 Subject: [PATCH 8/9] keep nmaxgr --- source/source_hamilt/module_surchem/cal_vcav.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index 4212a96816..85dded868f 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -40,7 +40,7 @@ void lapl_rho(const double& tpiba2, ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx); - std::complex *aux = new std::complex[rho_basis->nrxx]; + std::complex *aux = new std::complex[rho_basis->nmaxgr]; for (int ig = 0; ig < rho_basis->npw; ig++) { @@ -50,7 +50,7 @@ void lapl_rho(const double& tpiba2, for(int i = 0 ; i < 3 ; ++i) { - ModuleBase::GlobalFunc::ZEROS(aux, rho_basis->nrxx); + 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); From 0056c3625c752987710e8486f685506d68d69aed Mon Sep 17 00:00:00 2001 From: LKFEIYI <38131547+LKFEIYI@users.noreply.github.com> Date: Mon, 5 Jan 2026 19:26:20 -0700 Subject: [PATCH 9/9] roll back to original grad_rho --- .../source_hamilt/module_surchem/cal_vcav.cpp | 37 ++----------- .../source_hamilt/module_surchem/cal_vel.cpp | 54 +------------------ .../module_surchem/minimize_cg.cpp | 26 +-------- 3 files changed, 5 insertions(+), 112 deletions(-) diff --git a/source/source_hamilt/module_surchem/cal_vcav.cpp b/source/source_hamilt/module_surchem/cal_vcav.cpp index 85dded868f..b904e3ceea 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -4,33 +4,6 @@ #include "surchem.h" -static void local_grad_rho(const UnitCell& ucell, - const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, - ModuleBase::Vector3* grad_R, - std::complex* aux_G, - double* aux_R) -{ - - for (int i = 0; i < 3; ++i) - { - // 1. iG * rho(G) - for (int ig = 0; ig < rho_basis->npw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; - } - - // 2. FFT: G -> R - - rho_basis->recip2real(aux_G, aux_R); - - // 3. 2pi/a - for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_R[ir][i] = aux_R[ir] * ucell.tpiba; - } - } -} - - void lapl_rho(const double& tpiba2, const std::complex* rhog, double* lapn, @@ -111,10 +84,7 @@ void surchem::createcavity(const UnitCell& ucell, ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx); - std::complex *aux_G = new std::complex[rho_basis->npw]; - double *aux_R = new double[rho_basis->nrxx]; - - local_grad_rho(ucell, rho_basis, ps_totn, nablan, aux_G, aux_R); + XC_Functional::grad_rho(ps_totn, nablan, rho_basis, ucell.tpiba); // | \nabla n |^2 for (int ir = 0; ir < rho_basis->nrxx; ir++) @@ -164,7 +134,7 @@ void surchem::createcavity(const UnitCell& ucell, ModuleBase::Vector3 *ggn = new ModuleBase::Vector3[rho_basis->nrxx]; - local_grad_rho(ucell, rho_basis, inv_gn, ggn, aux_G, aux_R); + XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba); // -------------------------------------------------------- // add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav @@ -181,8 +151,7 @@ void surchem::createcavity(const UnitCell& ucell, } - delete[] aux_G; - delete[] aux_R; + delete[] nablan; delete[] nablan_2; delete[] sqrt_nablan_2; diff --git a/source/source_hamilt/module_surchem/cal_vel.cpp b/source/source_hamilt/module_surchem/cal_vel.cpp index 5022296b3b..c10ece46c3 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -4,29 +4,6 @@ #include "surchem.h" -static void local_grad_rho_vel(const UnitCell& ucell, - const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, - ModuleBase::Vector3* grad_R, - std::complex* aux_G, - double* aux_R) -{ - for (int i = 0; i < 3; ++i) - { - // 1. iG * rho(G) - for (int ig = 0; ig < rho_basis->npw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; - } - - // 2. FFT: G -> R - rho_basis->recip2real(aux_G, aux_R); - - // 3. 2pi/a - for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_R[ir][i] = aux_R[ir] * ucell.tpiba; - } - } -} void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis, double* eprime) { @@ -40,24 +17,6 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis } } -static void local_grad_rho_vel_tpiba(const double tpiba, - const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, - ModuleBase::Vector3* grad_R, - std::complex* aux_G, - double* aux_R) -{ - for (int i = 0; i < 3; ++i) - { - for (int ig = 0; ig < rho_basis->npw; ig++) { - aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; - } - rho_basis->recip2real(aux_G, aux_R); - for (int ir = 0; ir < rho_basis->nrxx; ir++) { - grad_R[ir][i] = aux_R[ir] * tpiba; - } - } -} void eps_pot(const double* PS_TOTN_real, const double& tpiba, @@ -80,17 +39,7 @@ void eps_pot(const double* PS_TOTN_real, double *phisq = new double[rho_basis->nrxx]; - - std::complex *aux_G = new std::complex[rho_basis->npw]; - double *aux_R = new double[rho_basis->nrxx]; - - - local_grad_rho_vel_tpiba(tpiba, rho_basis, phi, nabla_phi, aux_G, aux_R); - - - delete[] aux_G; - delete[] aux_R; - + XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba); for (int ir = 0; ir < rho_basis->nrxx; ir++) @@ -174,7 +123,6 @@ void surchem::cal_vel(const UnitCell& cell, this->Ael *= cell.omega / rho_basis->nxyz; // the 2nd item of tmp_Vel - // eps_pot internally calls shape_gradn and now the SAFE local_grad_rho_vel_tpiba eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot); for (int i = 0; i < rho_basis->nrxx; i++) diff --git a/source/source_hamilt/module_surchem/minimize_cg.cpp b/source/source_hamilt/module_surchem/minimize_cg.cpp index 3581f861e3..34489dadf9 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -177,30 +177,6 @@ void surchem::minimize_cg(const UnitCell& ucell, delete[] aux_grad_grad_phi_real; } -void helper_grad_rho(const UnitCell& ucell, - const ModulePW::PW_Basis* rho_basis, - const std::complex* rho_G, // G - ModuleBase::Vector3* grad_rho_R, // R - std::complex* buffer_G, - double* buffer_R) -{ - for (int i = 0; i < 3; ++i) - { - - for (int ig = 0; ig < rho_basis->npw; ig++) - { - buffer_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar[ig][i]; - } - - rho_basis->recip2real(buffer_G, buffer_R); - - - for (int ir = 0; ir < rho_basis->nrxx; ir++) - { - grad_rho_R[ir][i] = buffer_R[ir] * ucell.tpiba; - } - } -} void surchem::Leps2(const UnitCell& ucell, const ModulePW::PW_Basis* rho_basis, @@ -217,7 +193,7 @@ void surchem::Leps2(const UnitCell& ucell, double* aux_R) // size: nrxx { - helper_grad_rho(ucell, rho_basis, phi, grad_phi_R, aux_G, aux_R); + XC_Functional::grad_rho(phi, grad_phi_R, rho_basis, ucell.tpiba); for (int ir = 0; ir < rho_basis->nrxx; ir++)