From b087df9a4790b93372b8f05f34e4b852d469feda Mon Sep 17 00:00:00 2001 From: someone Date: Mon, 2 Mar 2026 18:46:13 +0800 Subject: [PATCH 1/9] INIT MGGA VELOCITY --- source/source_pw/module_pwdft/elecond.cpp | 21 ++++- source/source_pw/module_pwdft/op_pw_vel.cpp | 83 ++++++++++++++++++- source/source_pw/module_pwdft/op_pw_vel.h | 16 +++- .../source_pw/module_stodft/sto_elecond.cpp | 36 +++++++- 4 files changed, 147 insertions(+), 9 deletions(-) diff --git a/source/source_pw/module_pwdft/elecond.cpp b/source/source_pw/module_pwdft/elecond.cpp index a362a3088a9..3b1a82205cc 100644 --- a/source/source_pw/module_pwdft/elecond.cpp +++ b/source/source_pw/module_pwdft/elecond.cpp @@ -7,6 +7,8 @@ #include "source_base/parallel_device.h" #include "source_base/parallel_reduce.h" #include "source_estate/occupy.h" +#include "source_estate/module_pot/potential_new.h" +#include "source_base/module_device/types.h" #include "source_io/module_output/binstream.h" #include "source_io/module_parameter/parameter.h" @@ -86,7 +88,24 @@ void EleCond::KG(const int& smear_type, std::vector ct12(nt, 0); std::vector ct22(nt, 0); - hamilt::Velocity velop(this->p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, this->p_ucell, nonlocal); + using Real = typename GetTypeReal::type; + const Real* vtau_ptr = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->template get_vofk_smooth_data() + : nullptr; + const int vtau_col = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->get_vofk_smooth().nc + : 0; + const int vtau_row = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->get_vofk_smooth().nr + : 0; + hamilt::Velocity velop(this->p_wfcpw, + this->p_kv->isk.data(), + this->p_ppcell, + this->p_ucell, + nonlocal, + vtau_ptr, + vtau_col, + vtau_row); double decut = (wcut + fwhmin) / ModuleBase::Ry_to_eV; std::cout << "Recommended dt: " << 0.25 * M_PI / decut << " a.u." << std::endl; for (int ik = 0; ik < nk; ++ik) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index e0a8ed49ea5..fa8ab4606d4 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -3,6 +3,8 @@ #include "source_base/kernels/math_kernel_op.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" +#include "source_hamilt/module_xc/xc_functional.h" +#include "source_pw/module_pwdft/kernels/meta_op.h" namespace hamilt { @@ -11,7 +13,10 @@ Velocity::Velocity(const ModulePW::PW_Basis_K* wfcpw_in, const int* isk_in, pseudopot_cell_vnl* ppcell_in, const UnitCell* ucell_in, - const bool nonlocal_in) + const bool nonlocal_in, + const typename GetTypeReal::type* vtau_in, + const int vtau_col_in, + const int vtau_row_in) { if (wfcpw_in == nullptr || isk_in == nullptr || ppcell_in == nullptr || ucell_in == nullptr) { @@ -23,10 +28,16 @@ Velocity::Velocity(const ModulePW::PW_Basis_K* wfcpw_in, this->ucell = ucell_in; this->nonlocal = nonlocal_in; this->tpiba = ucell_in->tpiba; + this->vtau_ = vtau_in; + this->vtau_col_ = vtau_col_in; + this->vtau_row_ = vtau_row_in; if (this->nonlocal) { this->ppcell->initgradq_vnl(*this->ucell); } + // workspace for meta-GGA correction (size follows wfcpw grids) + resmem_complex_op()(porter1_, this->wfcpw->nmaxgr); + resmem_complex_op()(porter2_, this->wfcpw->nmaxgr); } template @@ -37,6 +48,8 @@ Velocity::~Velocity() delmem_var_op()(this->gz_); delmem_complex_op()(vkb_); delmem_complex_op()(gradvkb_); + delmem_complex_op()(porter1_); + delmem_complex_op()(porter2_); } template @@ -93,6 +106,7 @@ void Velocity::act(const psi::Psi, Device>* const int npw = this->wfcpw->npwk[this->ik]; const int max_npw = this->wfcpw->npwk_max; const int npol = psi_in->get_npol(); + using Real = typename GetTypeReal::type; std::vector gtmp_ptr = {this->gx_, this->gy_, this->gz_}; // ------------- @@ -110,6 +124,71 @@ void Velocity::act(const psi::Psi, Device>* } } + // --------------------------------------------- + // meta-GGA velocity correction: i[V_tau, r] + // V_tau implemented in ABACUS as -∇·(v_tau ∇), + // so [V_tau, r]_j = -[ ∂_j(v_tau ψ) + v_tau ∂_j ψ ]. + // The contribution to velocity is -i times the bracket. + // --------------------------------------------- + if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_func_type() == 3) + { + const int current_spin = (this->vtau_row_ > 1 && psi_in->get_npol() == 2) ? this->isk[this->ik] : 0; + const Real* vtau_spin = this->vtau_ + current_spin * this->vtau_col_; + Complex minus_i(0.0, -1.0); + for (int ib = 0; ib < n_npwx; ++ib) + { + const Complex* bandpsi = psi0 + ib * max_npw; + for (int id = 0; id < 3; ++id) + { + // term1: ∂_id (v_tau * psi) + this->wfcpw->recip_to_real(this->ctx, bandpsi, this->porter1_, this->ik); + ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter1_, this->porter1_, vtau_spin); + this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); + meta_pw_op()(this->ctx, + this->ik, + id, + npw, + max_npw, + this->tpiba, + this->wfcpw->template get_gcar_data(), + this->wfcpw->template get_kvec_c_data(), + this->porter1_, + this->porter2_, + false); + + // term2: v_tau * ∂_id psi (reuse porter1_) + meta_pw_op()(this->ctx, + this->ik, + id, + npw, + max_npw, + this->tpiba, + this->wfcpw->template get_gcar_data(), + this->wfcpw->template get_kvec_c_data(), + bandpsi, + this->porter1_, + false); + this->wfcpw->recip_to_real(this->ctx, this->porter1_, this->porter1_, this->ik); + ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter1_, this->porter1_, vtau_spin); + this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); + + // sum and apply -i + ModuleBase::vector_add_vector_op()(npw, + this->porter1_, + this->porter1_, + static_cast(1.0), + this->porter2_, + static_cast(1.0)); + ModuleBase::scal_op()(npw, &minus_i, this->porter1_, 1); + + // accumulate to vpsi (component id, band ib) + Complex* vpsi_slice = vpsi + id * n_npwx * max_npw + ib * max_npw; + Complex one = 1.0; + ModuleBase::axpy_op()(npw, &one, this->porter1_, 1, vpsi_slice, 1); + } + } + } + // --------------------------------------------- // i[V_NL, r] = (\nabla_q+\nabla_q')V_{NL}(q,q') // |\beta><\beta|\psi> @@ -334,4 +413,4 @@ template class Velocity; template class Velocity; #endif -} // namespace hamilt \ No newline at end of file +} // namespace hamilt diff --git a/source/source_pw/module_pwdft/op_pw_vel.h b/source/source_pw/module_pwdft/op_pw_vel.h index 211ce4007e0..b5aa6a42a45 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.h +++ b/source/source_pw/module_pwdft/op_pw_vel.h @@ -2,6 +2,7 @@ #define VELOCITY_PW_H #include "op_pw.h" #include "source_cell/unitcell.h" +#include "source_base/module_device/types.h" #include "source_pw/module_pwdft/vnl_pw.h" #include "source_basis/module_pw/pw_basis_k.h" namespace hamilt @@ -17,7 +18,10 @@ class Velocity const int* isk_in, pseudopot_cell_vnl* ppcell_in, const UnitCell* ucell_in, - const bool nonlocal_in = true + const bool nonlocal_in = true, + const typename GetTypeReal::type* vtau_in = nullptr, + const int vtau_col_in = 0, + const int vtau_row_in = 0 ); ~Velocity(); @@ -54,7 +58,13 @@ class Velocity int ik=0; double tpiba=0.0; - + const typename GetTypeReal::type* vtau_ = nullptr; ///< [CPU] meta-GGA vtau on real grid (nspin x nrxx_smooth) + int vtau_col_ = 0; ///< number of grid points per spin for vtau + int vtau_row_ = 0; ///< number of spin channels stored in vtau_ + std::complex* porter1_ = nullptr; ///< workspace on real grid + std::complex* porter2_ = nullptr; ///< workspace on real grid / recip grid + base_device::DEVICE_CPU* ctx = {}; + private: FPTYPE* gx_ = nullptr; ///<[Device, npwx] x component of G+K FPTYPE* gy_ = nullptr; ///<[Device, npwx] y component of G+K @@ -76,4 +86,4 @@ class Velocity using syncmem_complex_h2d_op = base_device::memory::synchronize_memory_op, Device, base_device::DEVICE_CPU>; }; } -#endif \ No newline at end of file +#endif diff --git a/source/source_pw/module_stodft/sto_elecond.cpp b/source/source_pw/module_stodft/sto_elecond.cpp index 4be5c41f266..6d6095f067f 100644 --- a/source/source_pw/module_stodft/sto_elecond.cpp +++ b/source/source_pw/module_stodft/sto_elecond.cpp @@ -8,6 +8,8 @@ #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/vector3.h" +#include "source_estate/module_pot/potential_new.h" +#include "source_base/module_device/types.h" #include "source_io/module_parameter/parameter.h" #include "sto_tool.h" @@ -615,8 +617,37 @@ void Sto_EleCond::sKG(const int& smear_type, // ik loop ModuleBase::timer::start("Sto_EleCond", "kloop"); - hamilt::Velocity velop(this->p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, this->p_ucell, nonlocal); - hamilt::Velocity low_velop(this->p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, this->p_ucell, nonlocal); + using Real = typename GetTypeReal::type; + using LowReal = typename GetTypeReal::type; + const Real* vtau_ptr = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->template get_vofk_smooth_data() + : nullptr; + const LowReal* vtau_ptr_low = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->template get_vofk_smooth_data() + : nullptr; + const int vtau_col = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->get_vofk_smooth().nc + : 0; + const int vtau_row = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->get_vofk_smooth().nr + : 0; + + hamilt::Velocity velop(this->p_wfcpw, + this->p_kv->isk.data(), + this->p_ppcell, + this->p_ucell, + nonlocal, + vtau_ptr, + vtau_col, + vtau_row); + hamilt::Velocity low_velop(this->p_wfcpw, + this->p_kv->isk.data(), + this->p_ppcell, + this->p_ucell, + nonlocal, + vtau_ptr_low, + vtau_col, + vtau_row); for (int ik = 0; ik < nk; ++ik) { velop.init(ik); @@ -1079,4 +1110,3 @@ template class Sto_EleCond; #if ((defined __CUDA) || (defined __ROCM)) template class Sto_EleCond; #endif - From a8556d582ec29f6bfdf50a1166199577656c10ca Mon Sep 17 00:00:00 2001 From: someone Date: Sat, 11 Apr 2026 17:38:17 +0800 Subject: [PATCH 2/9] Test --- source/source_pw/module_pwdft/op_pw_vel.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index fa8ab4606d4..4fc9bccc54a 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -123,7 +123,11 @@ void Velocity::act(const psi::Psi, Device>* tmpsi_in += max_npw; } } - + if (this->ppcell->nkb <= 0 || !this->nonlocal) + { + ModuleBase::timer::tick("Operator", "Velocity"); + return; + } // --------------------------------------------- // meta-GGA velocity correction: i[V_tau, r] // V_tau implemented in ABACUS as -∇·(v_tau ∇), @@ -188,7 +192,8 @@ void Velocity::act(const psi::Psi, Device>* } } } - + ModuleBase::timer::tick("Operator", "Velocity"); + return; // --------------------------------------------- // i[V_NL, r] = (\nabla_q+\nabla_q')V_{NL}(q,q') // |\beta><\beta|\psi> From ff4867e41718f86ed9eb6932ba496ebb03867ee0 Mon Sep 17 00:00:00 2001 From: someone Date: Sat, 11 Apr 2026 18:13:57 +0800 Subject: [PATCH 3/9] Revert "Test" This reverts commit bed1ff65afc50546b788691d65e50bfcb12e18b3. --- source/source_pw/module_pwdft/op_pw_vel.cpp | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index 4fc9bccc54a..fa8ab4606d4 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -123,11 +123,7 @@ void Velocity::act(const psi::Psi, Device>* tmpsi_in += max_npw; } } - if (this->ppcell->nkb <= 0 || !this->nonlocal) - { - ModuleBase::timer::tick("Operator", "Velocity"); - return; - } + // --------------------------------------------- // meta-GGA velocity correction: i[V_tau, r] // V_tau implemented in ABACUS as -∇·(v_tau ∇), @@ -192,8 +188,7 @@ void Velocity::act(const psi::Psi, Device>* } } } - ModuleBase::timer::tick("Operator", "Velocity"); - return; + // --------------------------------------------- // i[V_NL, r] = (\nabla_q+\nabla_q')V_{NL}(q,q') // |\beta><\beta|\psi> From 4fb585abbcd4c669d2d8185fbd1fd80d0ccf601c Mon Sep 17 00:00:00 2001 From: someone Date: Sat, 11 Apr 2026 23:27:17 +0800 Subject: [PATCH 4/9] Fix times 2 --- source/source_pw/module_pwdft/op_pw_vel.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index fa8ab4606d4..fcbb303ab5b 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -126,15 +126,18 @@ void Velocity::act(const psi::Psi, Device>* // --------------------------------------------- // meta-GGA velocity correction: i[V_tau, r] - // V_tau implemented in ABACUS as -∇·(v_tau ∇), - // so [V_tau, r]_j = -[ ∂_j(v_tau ψ) + v_tau ∂_j ψ ]. + // Libxc returns d(exc)/d(tau) with tau = |∇ψ|^2 / 2, while the PW path + // applies the operator through the unscaled |∇ψ|^2 channel. Restrict the + // fix to the velocity path by adding the missing 1/2 here: + // V_tau = -(1/2) ∇·(v_tau ∇), so + // [V_tau, r]_j = -(1/2) [ ∂_j(v_tau ψ) + v_tau ∂_j ψ ]. // The contribution to velocity is -i times the bracket. // --------------------------------------------- if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_func_type() == 3) { const int current_spin = (this->vtau_row_ > 1 && psi_in->get_npol() == 2) ? this->isk[this->ik] : 0; const Real* vtau_spin = this->vtau_ + current_spin * this->vtau_col_; - Complex minus_i(0.0, -1.0); + Complex minus_half_i(0.0, -0.5); for (int ib = 0; ib < n_npwx; ++ib) { const Complex* bandpsi = psi0 + ib * max_npw; @@ -172,14 +175,14 @@ void Velocity::act(const psi::Psi, Device>* ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter1_, this->porter1_, vtau_spin); this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); - // sum and apply -i + // sum and apply -(i/2) ModuleBase::vector_add_vector_op()(npw, this->porter1_, this->porter1_, static_cast(1.0), this->porter2_, static_cast(1.0)); - ModuleBase::scal_op()(npw, &minus_i, this->porter1_, 1); + ModuleBase::scal_op()(npw, &minus_half_i, this->porter1_, 1); // accumulate to vpsi (component id, band ib) Complex* vpsi_slice = vpsi + id * n_npwx * max_npw + ib * max_npw; From 13f16d6cb9488321a7a9f2ba6b926b5fb6c4a0f6 Mon Sep 17 00:00:00 2001 From: someone Date: Sun, 7 Jun 2026 14:18:51 +0800 Subject: [PATCH 5/9] Fix(pw): address meta-gga velocity review comments --- source/source_pw/module_pwdft/op_pw_vel.cpp | 66 ++++++++++++--------- source/source_pw/module_pwdft/op_pw_vel.h | 6 +- 2 files changed, 42 insertions(+), 30 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index fcbb303ab5b..8dd2c5306bc 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -35,9 +35,6 @@ Velocity::Velocity(const ModulePW::PW_Basis_K* wfcpw_in, { this->ppcell->initgradq_vnl(*this->ucell); } - // workspace for meta-GGA correction (size follows wfcpw grids) - resmem_complex_op()(porter1_, this->wfcpw->nmaxgr); - resmem_complex_op()(porter2_, this->wfcpw->nmaxgr); } template @@ -125,28 +122,48 @@ void Velocity::act(const psi::Psi, Device>* } // --------------------------------------------- - // meta-GGA velocity correction: i[V_tau, r] + // meta-GGA velocity correction // Libxc returns d(exc)/d(tau) with tau = |∇ψ|^2 / 2, while the PW path // applies the operator through the unscaled |∇ψ|^2 channel. Restrict the // fix to the velocity path by adding the missing 1/2 here: // V_tau = -(1/2) ∇·(v_tau ∇), so // [V_tau, r]_j = -(1/2) [ ∂_j(v_tau ψ) + v_tau ∂_j ψ ]. - // The contribution to velocity is -i times the bracket. + // With the velocity convention used here, the contribution added below is + // -i/2 times the bracketed sum. // --------------------------------------------- if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_func_type() == 3) { - const int current_spin = (this->vtau_row_ > 1 && psi_in->get_npol() == 2) ? this->isk[this->ik] : 0; + if (this->porter1_ == nullptr) + { + resmem_complex_op()(this->porter1_, this->wfcpw->nmaxgr); + } + if (this->porter2_ == nullptr) + { + resmem_complex_op()(this->porter2_, this->wfcpw->nmaxgr); + } + int current_spin = 0; + if (this->vtau_row_ > 1) + { + current_spin = this->isk[this->ik]; + if (current_spin < 0 || current_spin >= this->vtau_row_) + { + ModuleBase::WARNING_QUIT("Velocity", "invalid spin index for meta-GGA velocity correction"); + } + } const Real* vtau_spin = this->vtau_ + current_spin * this->vtau_col_; Complex minus_half_i(0.0, -0.5); for (int ib = 0; ib < n_npwx; ++ib) { const Complex* bandpsi = psi0 + ib * max_npw; + this->wfcpw->recip_to_real(this->ctx, bandpsi, this->porter1_, this->ik); + ModuleBase::vector_mul_vector_op()(this->vtau_col_, + this->porter1_, + this->porter1_, + vtau_spin); + this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); for (int id = 0; id < 3; ++id) { // term1: ∂_id (v_tau * psi) - this->wfcpw->recip_to_real(this->ctx, bandpsi, this->porter1_, this->ik); - ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter1_, this->porter1_, vtau_spin); - this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); meta_pw_op()(this->ctx, this->ik, id, @@ -158,8 +175,12 @@ void Velocity::act(const psi::Psi, Device>* this->porter1_, this->porter2_, false); + ModuleBase::scal_op()(npw, &minus_half_i, this->porter2_, 1); + Complex* vpsi_slice = vpsi + id * n_npwx * max_npw + ib * max_npw; + Complex one = 1.0; + ModuleBase::axpy_op()(npw, &one, this->porter2_, 1, vpsi_slice, 1); - // term2: v_tau * ∂_id psi (reuse porter1_) + // term2: v_tau * ∂_id psi meta_pw_op()(this->ctx, this->ik, id, @@ -169,25 +190,16 @@ void Velocity::act(const psi::Psi, Device>* this->wfcpw->template get_gcar_data(), this->wfcpw->template get_kvec_c_data(), bandpsi, - this->porter1_, + this->porter2_, false); - this->wfcpw->recip_to_real(this->ctx, this->porter1_, this->porter1_, this->ik); - ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter1_, this->porter1_, vtau_spin); - this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); - - // sum and apply -(i/2) - ModuleBase::vector_add_vector_op()(npw, - this->porter1_, - this->porter1_, - static_cast(1.0), + this->wfcpw->recip_to_real(this->ctx, this->porter2_, this->porter2_, this->ik); + ModuleBase::vector_mul_vector_op()(this->vtau_col_, this->porter2_, - static_cast(1.0)); - ModuleBase::scal_op()(npw, &minus_half_i, this->porter1_, 1); - - // accumulate to vpsi (component id, band ib) - Complex* vpsi_slice = vpsi + id * n_npwx * max_npw + ib * max_npw; - Complex one = 1.0; - ModuleBase::axpy_op()(npw, &one, this->porter1_, 1, vpsi_slice, 1); + this->porter2_, + vtau_spin); + this->wfcpw->real_to_recip(this->ctx, this->porter2_, this->porter2_, this->ik); + ModuleBase::scal_op()(npw, &minus_half_i, this->porter2_, 1); + ModuleBase::axpy_op()(npw, &one, this->porter2_, 1, vpsi_slice, 1); } } } diff --git a/source/source_pw/module_pwdft/op_pw_vel.h b/source/source_pw/module_pwdft/op_pw_vel.h index b5aa6a42a45..0aa31be0649 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.h +++ b/source/source_pw/module_pwdft/op_pw_vel.h @@ -61,9 +61,9 @@ class Velocity const typename GetTypeReal::type* vtau_ = nullptr; ///< [CPU] meta-GGA vtau on real grid (nspin x nrxx_smooth) int vtau_col_ = 0; ///< number of grid points per spin for vtau int vtau_row_ = 0; ///< number of spin channels stored in vtau_ - std::complex* porter1_ = nullptr; ///< workspace on real grid - std::complex* porter2_ = nullptr; ///< workspace on real grid / recip grid - base_device::DEVICE_CPU* ctx = {}; + mutable std::complex* porter1_ = nullptr; ///< workspace on real grid / recip grid + mutable std::complex* porter2_ = nullptr; ///< workspace on real grid / recip grid + Device* ctx = {}; private: FPTYPE* gx_ = nullptr; ///<[Device, npwx] x component of G+K From 036277ebb9db29efc09a7672b160567b59de9aba Mon Sep 17 00:00:00 2001 From: someone Date: Sun, 7 Jun 2026 14:32:02 +0800 Subject: [PATCH 6/9] Feature(pw): add meta-gga velocity conductivity switch --- docs/advanced/input_files/input-main.md | 8 ++++++++ .../source_io/module_ctrl/ctrl_output_pw.cpp | 1 + .../module_parameter/input_parameter.h | 1 + .../read_input_item_postprocess.cpp | 6 ++++++ source/source_io/test/read_input_ptest.cpp | 1 + source/source_io/test/support/INPUT | 5 +++++ source/source_pw/module_pwdft/elecond.cpp | 7 ++++--- source/source_pw/module_pwdft/elecond.h | 4 +++- source/source_pw/module_stodft/sto_elecond.cpp | 18 ++++++------------ 9 files changed, 35 insertions(+), 16 deletions(-) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index a91c6dd7530..a722057c512 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -486,6 +486,7 @@ - [cond\_smear](#cond_smear) - [cond\_fwhm](#cond_fwhm) - [cond\_nonlocal](#cond_nonlocal) + - [cond\_mgga\_vel](#cond_mgga_vel) - [Implicit solvation model](#implicit-solvation-model) - [imp\_sol](#imp_sol) - [eb\_k](#eb_k) @@ -4314,6 +4315,13 @@ - False: . - **Default**: True +### cond_mgga_vel + +- **Type**: Boolean +- **Availability**: [basis_type](#basis_type) = `pw` +- **Description**: Whether to include the meta-GGA velocity correction from the $v_\tau$ term when calculating velocity matrix $\bra{\psi_i}\hat{v}\ket{\psi_j}$. +- **Default**: True + [back to top](#full-list-of-input-keywords) ## Implicit solvation model diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index 580bc314efb..57e4cabf2e4 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -334,6 +334,7 @@ void ModuleIO::ctrl_runner_pw(UnitCell& ucell, inp.cond_dw, inp.cond_dt, inp.cond_nonlocal, + inp.cond_mgga_vel, pelec->wg); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index acbfeb62bf4..96a7cc57547 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -451,6 +451,7 @@ struct Input_para int cond_smear = 1; ///< smearing method for conductivities 1: Gaussian 2: Lorentzian double cond_fwhm = 0.4; ///< FWHM for conductivities bool cond_nonlocal = true; ///< if calculate nonlocal effects + bool cond_mgga_vel = true; ///< if calculate meta-GGA velocity correction bool berry_phase = false; ///< berry phase calculation: calculate berry phase or not int gdir = 3; ///< berry phase calculation: calculate the polarization in diff --git a/source/source_io/module_parameter/read_input_item_postprocess.cpp b/source/source_io/module_parameter/read_input_item_postprocess.cpp index b9a6e8f6ba1..aad6f47a85a 100644 --- a/source/source_io/module_parameter/read_input_item_postprocess.cpp +++ b/source/source_io/module_parameter/read_input_item_postprocess.cpp @@ -275,6 +275,12 @@ void ReadInput::item_postprocess() read_sync_bool(input.cond_nonlocal); this->add_item(item); } + { + Input_Item item("cond_mgga_vel"); + item.annotation = "Meta-GGA velocity correction for conductivities"; + read_sync_bool(input.cond_mgga_vel); + this->add_item(item); + } // berry_wannier { diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index 615757b112d..05f5b1e75f6 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -73,6 +73,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.cond_dtbatch, 2); EXPECT_DOUBLE_EQ(param.inp.cond_fwhm, 0.3); EXPECT_TRUE(param.inp.cond_nonlocal); + EXPECT_TRUE(param.inp.cond_mgga_vel); EXPECT_FALSE(param.inp.berry_phase); EXPECT_EQ(param.inp.ocp_kb.size(), 2); EXPECT_EQ(param.inp.ocp_kb[0], 1); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index 799c2e7a318..6915f667376 100644 --- a/source/source_io/test/support/INPUT +++ b/source/source_io/test/support/INPUT @@ -91,6 +91,7 @@ cond_dt 0.07 #control the t interval cond_dtbatch 2 #control dt batch cond_fwhm 0.3 #FWHM for conductivities cond_nonlocal 1 #Nonlocal effects for conductivities +cond_mgga_vel 1 #Meta-GGA velocity correction for conductivities #Parameters (4.Relaxation) ks_solver genelpa #cg; dav; lapack; genelpa; scalapack_gvx; cusolver @@ -391,3 +392,7 @@ nsc 50 #Maximal number of spin-constrained iteration nsc_min 4 #Minimum number of spin-constrained iteration alpha_trial 0.02 #Initial trial step size for lambda in eV/uB^2 sccut 4 #Maximal step size for lambda in eV/uB + +#Parameters (23. Time-dependent orbital-free DFT) +of_cd 0 #0: no CD potential; 1: add CD potential +of_mCD_alpha 1.0 # parameter of modified CD potential diff --git a/source/source_pw/module_pwdft/elecond.cpp b/source/source_pw/module_pwdft/elecond.cpp index 3b1a82205cc..80f5aa93362 100644 --- a/source/source_pw/module_pwdft/elecond.cpp +++ b/source/source_pw/module_pwdft/elecond.cpp @@ -53,6 +53,7 @@ void EleCond::KG(const int& smear_type, const double& dw_in, const double& dt_in, const bool& nonlocal, + const bool& mgga_vel, ModuleBase::matrix& wg) { //----------------------------------------------------------- @@ -89,13 +90,13 @@ void EleCond::KG(const int& smear_type, std::vector ct22(nt, 0); using Real = typename GetTypeReal::type; - const Real* vtau_ptr = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + const Real* vtau_ptr = (mgga_vel && this->p_elec != nullptr && this->p_elec->pot != nullptr) ? this->p_elec->pot->template get_vofk_smooth_data() : nullptr; - const int vtau_col = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + const int vtau_col = (mgga_vel && this->p_elec != nullptr && this->p_elec->pot != nullptr) ? this->p_elec->pot->get_vofk_smooth().nc : 0; - const int vtau_row = (this->p_elec != nullptr && this->p_elec->pot != nullptr) + const int vtau_row = (mgga_vel && this->p_elec != nullptr && this->p_elec->pot != nullptr) ? this->p_elec->pot->get_vofk_smooth().nr : 0; hamilt::Velocity velop(this->p_wfcpw, diff --git a/source/source_pw/module_pwdft/elecond.h b/source/source_pw/module_pwdft/elecond.h index 83a4a85d253..4e58a83f588 100644 --- a/source/source_pw/module_pwdft/elecond.h +++ b/source/source_pw/module_pwdft/elecond.h @@ -34,6 +34,7 @@ class EleCond * @param dw_in \omega step * @param dt_in time step * @param nonlocal whether to include the nonlocal potential corrections for velocity operator + * @param mgga_vel whether to include the meta-GGA velocity correction * @param wg wg(ik,ib) occupation for the ib-th band in the ik-th kpoint */ void KG(const int& smear_type, @@ -42,6 +43,7 @@ class EleCond const double& dw_in, const double& dt_in, const bool& nonlocal, + const bool& mgga_vel, ModuleBase::matrix& wg); protected: @@ -99,4 +101,4 @@ class EleCond double* ct22); }; -#endif // ELECOND_H \ No newline at end of file +#endif // ELECOND_H diff --git a/source/source_pw/module_stodft/sto_elecond.cpp b/source/source_pw/module_stodft/sto_elecond.cpp index 6d6095f067f..505ced5549b 100644 --- a/source/source_pw/module_stodft/sto_elecond.cpp +++ b/source/source_pw/module_stodft/sto_elecond.cpp @@ -619,18 +619,12 @@ void Sto_EleCond::sKG(const int& smear_type, ModuleBase::timer::start("Sto_EleCond", "kloop"); using Real = typename GetTypeReal::type; using LowReal = typename GetTypeReal::type; - const Real* vtau_ptr = (this->p_elec != nullptr && this->p_elec->pot != nullptr) - ? this->p_elec->pot->template get_vofk_smooth_data() - : nullptr; - const LowReal* vtau_ptr_low = (this->p_elec != nullptr && this->p_elec->pot != nullptr) - ? this->p_elec->pot->template get_vofk_smooth_data() - : nullptr; - const int vtau_col = (this->p_elec != nullptr && this->p_elec->pot != nullptr) - ? this->p_elec->pot->get_vofk_smooth().nc - : 0; - const int vtau_row = (this->p_elec != nullptr && this->p_elec->pot != nullptr) - ? this->p_elec->pot->get_vofk_smooth().nr - : 0; + // STO meta-GGA/SCAN is not implemented yet, so keep the meta-GGA velocity + // correction disabled for stochastic conductivity for now. + const Real* vtau_ptr = nullptr; + const LowReal* vtau_ptr_low = nullptr; + const int vtau_col = 0; + const int vtau_row = 0; hamilt::Velocity velop(this->p_wfcpw, this->p_kv->isk.data(), From c891000a34b0bceb6bc88a8827eda04390e9cf69 Mon Sep 17 00:00:00 2001 From: someone Date: Sun, 7 Jun 2026 14:38:07 +0800 Subject: [PATCH 7/9] Docs(pw): clarify meta-gga velocity correction --- source/source_pw/module_pwdft/op_pw_vel.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index 8dd2c5306bc..da2740c24bb 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -123,13 +123,13 @@ void Velocity::act(const psi::Psi, Device>* // --------------------------------------------- // meta-GGA velocity correction - // Libxc returns d(exc)/d(tau) with tau = |∇ψ|^2 / 2, while the PW path - // applies the operator through the unscaled |∇ψ|^2 channel. Restrict the - // fix to the velocity path by adding the missing 1/2 here: - // V_tau = -(1/2) ∇·(v_tau ∇), so - // [V_tau, r]_j = -(1/2) [ ∂_j(v_tau ψ) + v_tau ∂_j ψ ]. - // With the velocity convention used here, the contribution added below is - // -i/2 times the bracketed sum. + // V_tau = -(1/2) div(v_tau grad), whose plane-wave matrix element is + // = (1/2) v_tau(G-G') (k+G)·(k+G'). + // Therefore + // i[V_tau, r_alpha]_{G,G'} = + // (1/2) v_tau(G-G') [2k_alpha + G_alpha + G'_alpha]. + // In real space this is implemented as + // -i/2 [d_alpha(v_tau psi) + v_tau d_alpha psi]. // --------------------------------------------- if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_func_type() == 3) { From 266d82d3abad9aef6687c92409f78fb29ebfe5ac Mon Sep 17 00:00:00 2001 From: someone Date: Thu, 25 Jun 2026 23:17:22 +0800 Subject: [PATCH 8/9] Fix(pw): guard meta-gga velocity data --- source/source_pw/module_pwdft/elecond.cpp | 6 ++++++ source/source_pw/module_pwdft/op_pw_vel.cpp | 2 +- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/source/source_pw/module_pwdft/elecond.cpp b/source/source_pw/module_pwdft/elecond.cpp index 80f5aa93362..7e5bc4d4c75 100644 --- a/source/source_pw/module_pwdft/elecond.cpp +++ b/source/source_pw/module_pwdft/elecond.cpp @@ -8,6 +8,7 @@ #include "source_base/parallel_reduce.h" #include "source_estate/occupy.h" #include "source_estate/module_pot/potential_new.h" +#include "source_hamilt/module_xc/xc_functional.h" #include "source_base/module_device/types.h" #include "source_io/module_output/binstream.h" #include "source_io/module_parameter/parameter.h" @@ -99,6 +100,11 @@ void EleCond::KG(const int& smear_type, const int vtau_row = (mgga_vel && this->p_elec != nullptr && this->p_elec->pot != nullptr) ? this->p_elec->pot->get_vofk_smooth().nr : 0; + if (mgga_vel && XC_Functional::get_ked_flag() && (vtau_ptr == nullptr || vtau_col <= 0 || vtau_row <= 0)) + { + ModuleBase::WARNING_QUIT("EleCond::KG", + "meta-GGA velocity correction is requested, but v_tau data is unavailable"); + } hamilt::Velocity velop(this->p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index da2740c24bb..00c2f406b81 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -131,7 +131,7 @@ void Velocity::act(const psi::Psi, Device>* // In real space this is implemented as // -i/2 [d_alpha(v_tau psi) + v_tau d_alpha psi]. // --------------------------------------------- - if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_func_type() == 3) + if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_ked_flag()) { if (this->porter1_ == nullptr) { From e1145e93fe920f42a254db408aa4548090927da5 Mon Sep 17 00:00:00 2001 From: someone Date: Fri, 26 Jun 2026 00:16:57 +0800 Subject: [PATCH 9/9] Style(pw): use ASCII in velocity comments --- source/source_pw/module_pwdft/op_pw_vel.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/source/source_pw/module_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index 00c2f406b81..d2e6bb8109b 100644 --- a/source/source_pw/module_pwdft/op_pw_vel.cpp +++ b/source/source_pw/module_pwdft/op_pw_vel.cpp @@ -124,12 +124,12 @@ void Velocity::act(const psi::Psi, Device>* // --------------------------------------------- // meta-GGA velocity correction // V_tau = -(1/2) div(v_tau grad), whose plane-wave matrix element is - // = (1/2) v_tau(G-G') (k+G)·(k+G'). + // = (1/2) v_tau(G-G') (k+G) dot (k+G'). // Therefore - // i[V_tau, r_alpha]_{G,G'} = + // i[V_tau, r_\alpha]_{G,G'} = // (1/2) v_tau(G-G') [2k_alpha + G_alpha + G'_alpha]. // In real space this is implemented as - // -i/2 [d_alpha(v_tau psi) + v_tau d_alpha psi]. + // -i/2 [\partial_\alpha(v_tau psi) + v_tau \partial_\alpha psi]. // --------------------------------------------- if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_ked_flag()) { @@ -163,7 +163,7 @@ void Velocity::act(const psi::Psi, Device>* this->wfcpw->real_to_recip(this->ctx, this->porter1_, this->porter1_, this->ik); for (int id = 0; id < 3; ++id) { - // term1: ∂_id (v_tau * psi) + // term1: partial_id (v_tau * psi) meta_pw_op()(this->ctx, this->ik, id, @@ -180,7 +180,7 @@ void Velocity::act(const psi::Psi, Device>* Complex one = 1.0; ModuleBase::axpy_op()(npw, &one, this->porter2_, 1, vpsi_slice, 1); - // term2: v_tau * ∂_id psi + // term2: v_tau * partial_id psi meta_pw_op()(this->ctx, this->ik, id,