support force metagga calculation

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp

@@ -124,30 +124,27 @@ void Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>::contributeH
     ModuleBase::timer::tick("Veff", "contributeHR");
 
     std::vector<const double*> vr_eff(4, nullptr);
+    std::vector<const double*> vofk_eff(4, nullptr);
     for (int is = 0; is < 4; is++)
     {
-        const double* vr_eff1 = this->pot->get_effective_v(is);
+        vr_eff[is] = this->pot->get_effective_v(is);
         if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5)
         {
-            const double* vofk_eff1 = this->pot->get_effective_vofk(is);
-            Gint_inout inout(vr_eff1, vofk_eff1, is, Gint_Tools::job_type::vlocal_meta);
-            this->GK->cal_gint(&inout);
+            vofk_eff[is] = this->pot->get_effective_vofk(is);
+            if(is == 3)
+            {
+                ModuleGint::cal_gint_vl_metagga(vr_eff, vofk_eff, this->hR);
+            }
         }
         else
         {
-            vr_eff[is] = vr_eff1;
             if(is == 3)
             {
                 ModuleGint::cal_gint_vl(vr_eff, this->hR);
             }
         }
     }
 
-    if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5)
-    {
-        this->GK->transfer_pvpR(this->hR,this->ucell,this->gd);
-    }
-
     ModuleBase::timer::tick("Veff", "contributeHR");
     return;
 }

source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp

@@ -18,22 +18,20 @@ namespace PulayForceStress
         const bool& isstress,
         const bool& set_dmr_gint)
     {
-        if (set_dmr_gint) { gint.transfer_DM2DtoGrid(dm.get_DMR_vector()); }    // 2d block to grid
         const int nspin = PARAM.inp.nspin;
+        std::vector<const double*> vr_eff(nspin, nullptr);
+        std::vector<const double*> vofk_eff(nspin, nullptr);
         if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
         {
             for (int is = 0; is < nspin; ++is)
             {
-                const double* vr_eff1 = pot->get_effective_v(is);
-                const double* vofk_eff1 = nullptr;
-                vofk_eff1 = pot->get_effective_vofk(is);
-                Gint_inout inout(is, vr_eff1, vofk_eff1, isforce, isstress, &f, &s, Gint_Tools::job_type::force_meta);
-                gint.cal_gint(&inout);
+                vr_eff[is] = pot->get_effective_v(is);
+                vofk_eff[is] = pot->get_effective_vofk(is);
             }
+            ModuleGint::cal_gint_fvl_meta(nspin, vr_eff, vofk_eff, dm.get_DMR_vector(), isforce, isstress, &f, &s);
         }
         else
         {
-            std::vector<const double*> vr_eff(nspin, nullptr);
             for(int is = 0; is < nspin; ++is)
             {
                 vr_eff[is] = pot->get_effective_v(is);

source/module_hamilt_lcao/module_gint/CMakeLists.txt

@@ -39,9 +39,11 @@ list(APPEND objects
     new_grid_tech/gint_vl.cpp
     new_grid_tech/gint_vl_metagga.cpp
     new_grid_tech/gint_vl_nspin4.cpp
+    new_grid_tech/gint_vl_metagga_nspin4.cpp
     new_grid_tech/gint_rho.cpp
     new_grid_tech/gint_tau.cpp
     new_grid_tech/gint_fvl.cpp
+    new_grid_tech/gint_fvl_meta.cpp
     new_grid_tech/localcell_info.cpp
     new_grid_tech/phi_operator.cpp
     new_grid_tech/set_ddphi.cpp

source/module_hamilt_lcao/module_gint/new_grid_tech/gint_fvl_meta.cpp

@@ -0,0 +1,122 @@
+#include "module_base/array_pool.h"
+#include "module_base/global_function.h"
+#include "gint_fvl_meta.h"
+#include "gint_common.h"
+#include "phi_operator.h"
+
+namespace ModuleGint
+{
+
+void Gint_fvl_meta::cal_gint()
+{
+    init_DMRGint_();
+    transfer_DM_to_DMGint(gint_info_, DMR_vec_, DMRGint_vec_);
+    cal_fvl_svl_();
+}
+
+void Gint_fvl_meta::init_DMRGint_()
+{
+    DMRGint_vec_.resize(nspin_);
+    for (int is = 0; is < nspin_; is++)
+    {
+        DMRGint_vec_[is] = gint_info_->get_hr<double>();
+    }
+}
+
+void Gint_fvl_meta::cal_fvl_svl_()
+{
+#pragma omp parallel
+    {
+        PhiOperator phi_op;
+        ModuleBase::matrix* fvl_thread = nullptr;
+        ModuleBase::matrix* svl_thread = nullptr;
+        if(isforce_)
+        {
+            fvl_thread = new ModuleBase::matrix(*fvl_);
+            fvl_thread->zero_out();
+        }
+        if(isstress_)
+        {
+            svl_thread = new ModuleBase::matrix(*svl_);
+            svl_thread->zero_out();
+        }
+#pragma omp for schedule(dynamic)
+        for(const auto& biggrid: gint_info_->get_biggrids())
+        {
+            if(biggrid->get_atoms().size() == 0)
+            {
+                continue;
+            }
+            phi_op.set_bgrid(biggrid);
+            ModuleBase::Array_Pool<double> phi(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_x(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_y(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_z(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_xx(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_xy(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_xz(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_yy(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_yz(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> ddphi_zz(phi_op.get_rows(), phi_op.get_cols());
+            phi_op.set_phi_dphi(phi.get_ptr_1D(), dphi_x.get_ptr_1D(), dphi_y.get_ptr_1D(), dphi_z.get_ptr_1D());
+            phi_op.set_ddphi(ddphi_xx.get_ptr_1D(), ddphi_xy.get_ptr_1D(), ddphi_xz.get_ptr_1D(),
+                             ddphi_yy.get_ptr_1D(), ddphi_yz.get_ptr_1D(), ddphi_zz.get_ptr_1D());
+            ModuleBase::Array_Pool<double> phi_vldr3(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> phi_vldr3_DM(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_x_vldr3(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_y_vldr3(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_z_vldr3(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_x_vldr3_DM(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_y_vldr3_DM(phi_op.get_rows(), phi_op.get_cols());
+            ModuleBase::Array_Pool<double> dphi_z_vldr3_DM(phi_op.get_rows(), phi_op.get_cols());
+            for (int is = 0; is < nspin_; is++)
+            {
+                ModuleBase::zeros(phi_vldr3.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(phi_vldr3_DM.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_x_vldr3.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_y_vldr3.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_z_vldr3.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_x_vldr3_DM.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_y_vldr3_DM.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                ModuleBase::zeros(dphi_z_vldr3_DM.get_ptr_1D(), phi_op.get_rows()*phi_op.get_cols());
+                phi_op.phi_mul_vldr3(vr_eff_[is], dr3_, phi.get_ptr_2D(), phi_vldr3.get_ptr_2D());
+                phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_x.get_ptr_2D(), dphi_x_vldr3.get_ptr_2D());
+                phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_y.get_ptr_2D(), dphi_y_vldr3.get_ptr_2D());
+                phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_z.get_ptr_2D(), dphi_z_vldr3.get_ptr_2D());
+                phi_op.phi_mul_dm(phi_vldr3.get_ptr_2D(), *DMRGint_vec_[is], false, phi_vldr3_DM.get_ptr_2D());
+                phi_op.phi_mul_dm(dphi_x_vldr3.get_ptr_2D(), *DMRGint_vec_[is], false, dphi_x_vldr3_DM.get_ptr_2D());
+                phi_op.phi_mul_dm(dphi_y_vldr3.get_ptr_2D(), *DMRGint_vec_[is], false, dphi_y_vldr3_DM.get_ptr_2D());
+                phi_op.phi_mul_dm(dphi_z_vldr3.get_ptr_2D(), *DMRGint_vec_[is], false, dphi_z_vldr3_DM.get_ptr_2D());
+                if(isforce_)
+                {
+                    phi_op.phi_dot_dphi(phi_vldr3_DM.get_ptr_2D(), dphi_x.get_ptr_2D(), dphi_y.get_ptr_2D(), dphi_z.get_ptr_2D(), fvl_thread);
+                    phi_op.phi_dot_dphi(dphi_x_vldr3_DM.get_ptr_2D(), ddphi_xx.get_ptr_2D(), ddphi_xy.get_ptr_2D(), ddphi_xz.get_ptr_2D(), fvl_thread);
+                    phi_op.phi_dot_dphi(dphi_y_vldr3_DM.get_ptr_2D(), ddphi_xy.get_ptr_2D(), ddphi_yy.get_ptr_2D(), ddphi_yz.get_ptr_2D(), fvl_thread);
+                    phi_op.phi_dot_dphi(dphi_z_vldr3_DM.get_ptr_2D(), ddphi_xz.get_ptr_2D(), ddphi_yz.get_ptr_2D(), ddphi_zz.get_ptr_2D(), fvl_thread);
+                }
+                if(isstress_)
+                {
+                    phi_op.phi_dot_dphi_r(phi_vldr3_DM.get_ptr_2D(), dphi_x.get_ptr_2D(), dphi_y.get_ptr_2D(), dphi_z.get_ptr_2D(), svl_thread);
+                    phi_op.phi_dot_dphi_r(dphi_x_vldr3_DM.get_ptr_2D(), ddphi_xx.get_ptr_2D(), ddphi_xy.get_ptr_2D(), ddphi_xz.get_ptr_2D(), svl_thread);
+                    phi_op.phi_dot_dphi_r(dphi_y_vldr3_DM.get_ptr_2D(), ddphi_xy.get_ptr_2D(), ddphi_yy.get_ptr_2D(), ddphi_yz.get_ptr_2D(), svl_thread);
+                    phi_op.phi_dot_dphi_r(dphi_z_vldr3_DM.get_ptr_2D(), ddphi_xz.get_ptr_2D(), ddphi_yz.get_ptr_2D(), ddphi_zz.get_ptr_2D(), svl_thread);
+                }
+            }
+        }
+#pragma omp critical
+        {
+            if(isforce_)
+            {
+                fvl_[0] += fvl_thread[0];
+                delete fvl_thread;
+            }
+            if(isstress_)
+            {
+                svl_[0] += svl_thread[0];
+                delete svl_thread;
+            }
+        }
+    }
+}
+
+} // namespace ModuleGint

source/module_hamilt_lcao/module_gint/new_grid_tech/gint_fvl_meta.h

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <memory>
+#include <vector>
+#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
+#include "module_base/matrix.h"
+#include "gint.h"
+#include "gint_info.h"
+
+namespace ModuleGint
+{
+class Gint_fvl_meta : public Gint
+{
+    public:
+    Gint_fvl_meta(
+        const int nspin,
+        const std::vector<const double*>& vr_eff,
+        const std::vector<const double*>& vofk,
+        const std::vector<HContainer<double>*>& DMR_vec,
+        const bool isforce,
+        const bool isstress,
+        ModuleBase::matrix* fvl,
+        ModuleBase::matrix* svl)
+        : nspin_(nspin), vr_eff_(vr_eff), vofk_(vofk), DMR_vec_(DMR_vec),
+          isforce_(isforce), isstress_(isstress), fvl_(fvl), svl_(svl),
+          dr3_(gint_info_->get_mgrid_volume()) {};
+
+    void cal_gint() override;
+
+    private:
+    void init_DMRGint_();
+
+    void cal_fvl_svl_();
+
+    // input
+    const int nspin_;
+    std::vector<const double*> vr_eff_;
+    std::vector<const double*> vofk_;
+    std::vector<HContainer<double>*> DMR_vec_;
+    const bool isforce_;
+    const bool isstress_;
+
+    // output
+    ModuleBase::matrix* fvl_;
+    ModuleBase::matrix* svl_;
+
+    // intermediate variables
+    std::vector<std::shared_ptr<HContainer<double>>> DMRGint_vec_;
+
+    double dr3_;
+};
+
+} // namespace ModuleGint

source/module_hamilt_lcao/module_gint/new_grid_tech/gint_interface.cpp

@@ -2,7 +2,9 @@
 #include "gint_vl.h"
 #include "gint_vl_metagga.h"
 #include "gint_vl_nspin4.h"
+#include "gint_vl_metagga_nspin4.h"
 #include "gint_fvl.h"
+#include "gint_fvl_meta.h"
 #include "gint_rho.h"
 #include "gint_tau.h"
 
@@ -34,6 +36,15 @@ void cal_gint_vl_metagga(
     gint_vl_metagga.cal_gint();
 }
 
+void cal_gint_vl_metagga(
+    std::vector<const double*> vr_eff,
+    std::vector<const double*> vofk,
+    HContainer<std::complex<double>>* hR)
+{
+    Gint_vl_metagga_nspin4 gint_vl_metagga_nspin4(vr_eff, vofk, hR);
+    gint_vl_metagga_nspin4.cal_gint();
+}
+
 void cal_gint_rho(
     const std::vector<HContainer<double>*>& DMR_vec,
     const int nspin,
@@ -65,4 +76,18 @@ void cal_gint_fvl(
     gint_fvl.cal_gint();
 }
 
+void cal_gint_fvl_meta(
+    const int nspin,
+    const std::vector<const double*>& vr_eff,
+    const std::vector<const double*>& vofk,
+    const std::vector<HContainer<double>*>& DMR_vec,
+    const bool isforce,
+    const bool isstress,
+    ModuleBase::matrix* fvl,
+    ModuleBase::matrix* svl)
+{
+    Gint_fvl_meta gint_fvl_meta(nspin, vr_eff, vofk, DMR_vec, isforce, isstress, fvl, svl);
+    gint_fvl_meta.cal_gint();
+}
+
 } // namespace ModuleGint

source/module_hamilt_lcao/module_gint/new_grid_tech/gint_interface.h

@@ -20,6 +20,11 @@ void cal_gint_vl_metagga(
     const double* vfork,
     HContainer<double>* hR);
 
+void cal_gint_vl_metagga(
+    std::vector<const double*> vr_eff,
+    std::vector<const double*> vofk,
+    HContainer<std::complex<double>>* hR);
+
 void cal_gint_rho(
     const std::vector<HContainer<double>*>& DMR_vec,
     const int nspin,
@@ -39,6 +44,16 @@ void cal_gint_fvl(
     ModuleBase::matrix* fvl,
     ModuleBase::matrix* svl);
 
+void cal_gint_fvl_meta(
+    const int nspin,
+    const std::vector<const double*>& vr_eff,
+    const std::vector<const double*>& vofk,
+    const std::vector<HContainer<double>*>& DMR_vec,
+    const bool isforce,
+    const bool isstress,
+    ModuleBase::matrix* fvl,
+    ModuleBase::matrix* svl);
+
 
 
 } // namespace ModuleGint