Skip to content

Fix memory leak in efield.cpp and optimize surchem module performance #6832

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
LKFEIYI wants to merge 10 commits into from
Closed

Fix memory leak in efield.cpp and optimize surchem module performance #6832

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
f474d9f
sol_refine
LKFEIYI Jan 4, 2026
474c403
fix void
LKFEIYI Jan 4, 2026
4d60da8
mem_leak_try
LKFEIYI Jan 5, 2026
896e63d
mem_leak
LKFEIYI Jan 5, 2026
3bd74b0
clear
LKFEIYI Jan 5, 2026
42b6949
Merge branch 'deepmodeling:develop' into develop
LKFEIYI Jan 5, 2026
dbd1a2b
update test
LKFEIYI Jan 5, 2026
bbec623
fix test
LKFEIYI Jan 5, 2026
4725616
keep nmaxgr
LKFEIYI Jan 5, 2026
0056c36
roll back to original grad_rho
LKFEIYI Jan 6, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion source/source_estate/module_pot/efield.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}

Expand Down
18 changes: 10 additions & 8 deletions source/source_estate/module_pot/pot_surchem.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<ModulePW::PW_Basis*>(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<ModulePW::PW_Basis*>(this->rho_basis_),
v_eff.nr,
chg->rho,
this->vlocal,
this->structure_factors_,
v_sol_correction);
v_eff += v_sol_correction;
}

private:
Expand Down
28 changes: 17 additions & 11 deletions source/source_hamilt/module_surchem/H_correction_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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;
Expand All @@ -58,5 +64,5 @@ ModuleBase::matrix surchem::v_correction(const UnitCell& cell,
delete[] total_n;

ModuleBase::timer::tick("surchem", "v_cor");
return v;
return;
}
2 changes: 1 addition & 1 deletion source/source_hamilt/module_surchem/cal_pseudo.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ void surchem::gauss_charge(const UnitCell& cell,
std::complex<double>* 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++)
Expand Down
106 changes: 55 additions & 51 deletions source/source_hamilt/module_surchem/cal_vcav.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "source_io/module_parameter/parameter.h"
#include "surchem.h"


void lapl_rho(const double& tpiba2,
const std::complex<double>* rhog,
double* lapn,
Expand All @@ -11,37 +12,41 @@ void lapl_rho(const double& tpiba2,
std::complex<double> *gdrtmpg = new std::complex<double>[rho_basis->npw];
ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);


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++) {
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
// exp(-(log(n/n_c))^2 /(2 sigma^2)) /(sigma * sqrt(2*pi) )/n
// ---------------------------------------------------------
// Helper 3: Shape function
// ---------------------------------------------------------
void shape_gradn(const std::complex<double>* 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);
Expand All @@ -59,13 +64,17 @@ void shape_gradn(const std::complex<double>* 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<double>* ps_totn,
double* vwork)
{
ModuleBase::Vector3<double> *nablan = new ModuleBase::Vector3<double>[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];
Expand All @@ -74,22 +83,18 @@ 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
// | \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;
Expand All @@ -100,48 +105,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<double> *inv_gn = new std::complex<double>[rho_basis->npw];
rho_basis->real2recip(sqrt_nablan_2, inv_gn);

// \nabla(1 / |\nabla n|), ggn in real space

ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];


XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba);

//-------------------------------------------------------------
// 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;
Expand All @@ -153,6 +150,8 @@ void surchem::createcavity(const UnitCell& ucell,
vwork[ir] = vwork[ir] * term1[ir] * PARAM.inp.tau;
}



delete[] nablan;
delete[] nablan_2;
delete[] sqrt_nablan_2;
Expand All @@ -162,10 +161,14 @@ void surchem::createcavity(const UnitCell& ucell,
delete[] ggn;
}

ModuleBase::matrix surchem::cal_vcav(const UnitCell& ucell,
const ModulePW::PW_Basis* rho_basis,
std::complex<double>* ps_totn,
int nspin)
// ---------------------------------------------------------
// Main Function: cal_vcav
// ---------------------------------------------------------
void surchem::cal_vcav(const UnitCell& ucell,
const ModulePW::PW_Basis* rho_basis,
std::complex<double>* ps_totn,
int nspin,
ModuleBase::matrix& v)
{
ModuleBase::TITLE("surchem", "cal_vcav");
ModuleBase::timer::tick("surchem", "cal_vcav");
Expand All @@ -175,12 +178,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
Expand All @@ -189,12 +192,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;
}
24 changes: 15 additions & 9 deletions source/source_hamilt/module_surchem/cal_vel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,10 @@
#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)
{

double epr_c = 1.0 / sqrt(ModuleBase::TWO_PI) / PARAM.inp.sigma_k;
double epr_z = 0;
double min = 1e-10;
Expand All @@ -16,6 +17,7 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
}
}


void eps_pot(const double* PS_TOTN_real,
const double& tpiba,
const std::complex<double>* phi,
Expand All @@ -36,9 +38,10 @@ void eps_pot(const double* PS_TOTN_real,
ModuleBase::Vector3<double> *nabla_phi = new ModuleBase::Vector3<double>[rho_basis->nrxx];
double *phisq = new double[rho_basis->nrxx];

// nabla phi

XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba);


for (int ir = 0; ir < rho_basis->nrxx; ir++)
{
phisq[ir] = pow(nabla_phi[ir].x, 2) + pow(nabla_phi[ir].y, 2) + pow(nabla_phi[ir].z, 2);
Expand All @@ -54,11 +57,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<double>* TOTN,
std::complex<double>* PS_TOTN,
int nspin)
void surchem::cal_vel(const UnitCell& cell,
const ModulePW::PW_Basis* rho_basis,
std::complex<double>* TOTN,
std::complex<double>* PS_TOTN,
int nspin,
ModuleBase::matrix& v)
{
ModuleBase::TITLE("surchem", "cal_vel");
ModuleBase::timer::tick("surchem", "cal_vel");
Expand Down Expand Up @@ -134,6 +138,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
Expand All @@ -143,6 +148,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);
}
}
}
Expand All @@ -158,5 +164,5 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
delete[] phi_tilda_R0;

ModuleBase::timer::tick("surchem", "cal_vel");
return Vel;
}
return;
}
Loading
Loading