diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e9dd89d7790..e9be4a71374 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) @@ -4322,6 +4323,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 34a548b576f..d554c44eb0a 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -452,6 +452,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 a362a3088a9..7e5bc4d4c75 100644 --- a/source/source_pw/module_pwdft/elecond.cpp +++ b/source/source_pw/module_pwdft/elecond.cpp @@ -7,6 +7,9 @@ #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_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" @@ -51,6 +54,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) { //----------------------------------------------------------- @@ -86,7 +90,29 @@ 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 = (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 = (mgga_vel && this->p_elec != nullptr && this->p_elec->pot != nullptr) + ? this->p_elec->pot->get_vofk_smooth().nc + : 0; + 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, + 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/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_pwdft/op_pw_vel.cpp b/source/source_pw/module_pwdft/op_pw_vel.cpp index e0a8ed49ea5..d2e6bb8109b 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,6 +28,9 @@ 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); @@ -37,6 +45,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 +103,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 +121,89 @@ 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) dot (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 [\partial_\alpha(v_tau psi) + v_tau \partial_\alpha psi]. + // --------------------------------------------- + if (this->vtau_ != nullptr && this->vtau_col_ > 0 && XC_Functional::get_ked_flag()) + { + 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: partial_id (v_tau * psi) + 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); + 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 * partial_id psi + 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->porter2_, + false); + this->wfcpw->recip_to_real(this->ctx, this->porter2_, this->porter2_, this->ik); + ModuleBase::vector_mul_vector_op()(this->vtau_col_, + this->porter2_, + 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); + } + } + } + // --------------------------------------------- // i[V_NL, r] = (\nabla_q+\nabla_q')V_{NL}(q,q') // |\beta><\beta|\psi> @@ -334,4 +428,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..0aa31be0649 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_ + 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 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..505ced5549b 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,31 @@ 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; + // 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(), + 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 +1104,3 @@ template class Sto_EleCond; #if ((defined __CUDA) || (defined __ROCM)) template class Sto_EleCond; #endif -