Feature: velocity matrix correction for meta-gga xc functionals (#7391)

* INIT MGGA VELOCITY

* Test

* Revert "Test"

This reverts commit bed1ff6.

* Fix times 2

* Fix(pw): address meta-gga velocity review comments

* Feature(pw): add meta-gga velocity conductivity switch

* Docs(pw): clarify meta-gga velocity correction

* Fix(pw): guard meta-gga velocity data

* Style(pw): use ASCII in velocity comments

Author: Flying-dragon-boxing

docs/advanced/input_files/input-main.md

@@ -502,6 +502,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)
@@ -4471,6 +4472,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

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);
     }
 

source/source_io/module_parameter/input_parameter.h

@@ -467,6 +467,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

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
     {

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);

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

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<FPTYPE, Device>::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<FPTYPE, Device>::KG(const int& smear_type,
     std::vector<double> ct12(nt, 0);
     std::vector<double> ct22(nt, 0);
 
-    hamilt::Velocity<FPTYPE, Device> velop(this->p_wfcpw, this->p_kv->isk.data(), this->p_ppcell, this->p_ucell, nonlocal);
+    using Real = typename GetTypeReal<FPTYPE>::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<Real>()
+                               : 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<FPTYPE, Device> 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)

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
+#endif // ELECOND_H

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<FPTYPE, Device>::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<FPTYPE>::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<FPTYPE, Device>::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<FPTYPE, Device>::~Velocity()
     delmem_var_op()(this->gz_);
     delmem_complex_op()(vkb_);
     delmem_complex_op()(gradvkb_);
+    delmem_complex_op()(porter1_);
+    delmem_complex_op()(porter2_);
 }
 
 template <typename FPTYPE, typename Device>
@@ -93,6 +103,7 @@ void Velocity<FPTYPE, Device>::act(const psi::Psi<std::complex<FPTYPE>, 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<FPTYPE>::type;
 
     std::vector<FPTYPE*> gtmp_ptr = {this->gx_, this->gy_, this->gz_};
     // -------------
@@ -110,6 +121,89 @@ void Velocity<FPTYPE, Device>::act(const psi::Psi<std::complex<FPTYPE>, Device>*
         }
     }
 
+    // ---------------------------------------------
+    // meta-GGA velocity correction
+    // V_tau = -(1/2) div(v_tau grad), whose plane-wave matrix element is
+    // <k+G|V_tau|k+G'> = (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<Complex, Device>()(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<Real, Device>()(this->ctx,
+                                           this->ik,
+                                           id,
+                                           npw,
+                                           max_npw,
+                                           this->tpiba,
+                                           this->wfcpw->template get_gcar_data<Real>(),
+                                           this->wfcpw->template get_kvec_c_data<Real>(),
+                                           this->porter1_,
+                                           this->porter2_,
+                                           false);
+                ModuleBase::scal_op<Real, Device>()(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<Complex, Device>()(npw, &one, this->porter2_, 1, vpsi_slice, 1);
+
+                // term2: v_tau * partial_id psi
+                meta_pw_op<Real, Device>()(this->ctx,
+                                           this->ik,
+                                           id,
+                                           npw,
+                                           max_npw,
+                                           this->tpiba,
+                                           this->wfcpw->template get_gcar_data<Real>(),
+                                           this->wfcpw->template get_kvec_c_data<Real>(),
+                                           bandpsi,
+                                           this->porter2_,
+                                           false);
+                this->wfcpw->recip_to_real(this->ctx, this->porter2_, this->porter2_, this->ik);
+                ModuleBase::vector_mul_vector_op<Complex, Device>()(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<Real, Device>()(npw, &minus_half_i, this->porter2_, 1);
+                ModuleBase::axpy_op<Complex, Device>()(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<double, base_device::DEVICE_GPU>;
 template class Velocity<float, base_device::DEVICE_GPU>;
 #endif
 
-} // namespace hamilt
+} // namespace hamilt

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<FPTYPE>::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<FPTYPE>::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<FPTYPE>* porter1_ = nullptr; ///< workspace on real grid / recip grid
+    mutable std::complex<FPTYPE>* 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<std::complex<FPTYPE>, Device, base_device::DEVICE_CPU>;
 };
 }
-#endif
+#endif