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_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..b904e3ceea 100644 --- a/source/source_hamilt/module_surchem/cal_vcav.cpp +++ b/source/source_hamilt/module_surchem/cal_vcav.cpp @@ -3,6 +3,7 @@ #include "source_io/module_parameter/parameter.h" #include "surchem.h" + void lapl_rho(const double& tpiba2, const std::complex* rhog, double* lapn, @@ -11,37 +12,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 + 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* 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 +64,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 +74,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 +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; @@ -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 *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); - //------------------------------------------------------------- - // 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 +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; @@ -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* ps_totn, - int nspin) +// --------------------------------------------------------- +// Main Function: cal_vcav +// --------------------------------------------------------- +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 +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 @@ -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; +} \ 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..c10ece46c3 100644 --- a/source/source_hamilt/module_surchem/cal_vel.cpp +++ b/source/source_hamilt/module_surchem/cal_vel.cpp @@ -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; @@ -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* phi, @@ -36,9 +38,10 @@ 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); + 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); @@ -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* TOTN, - std::complex* PS_TOTN, - int nspin) +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 +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 @@ -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); } } } @@ -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; +} \ 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..34489dadf9 100644 --- a/source/source_hamilt/module_surchem/minimize_cg.cpp +++ b/source/source_hamilt/module_surchem/minimize_cg.cpp @@ -8,16 +8,16 @@ 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 - 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 +29,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 +51,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 +73,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 +95,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 +111,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 +158,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 +167,84 @@ 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, + 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..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() 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, 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 65fe839aa2..a864f0d152 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); - - 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); + 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(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;