diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d6e74f898..d1487af838 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -236,6 +236,9 @@ if(ENABLE_LCAO) add_compile_definitions(__PEXSI) set(CMAKE_CXX_STANDARD 14) endif() + if(NEW_GINT) + add_compile_definitions(__NEW_GINT) + endif() else() set(ENABLE_DEEPKS OFF) set(ENABLE_LIBRI OFF) diff --git a/source/module_base/array_pool.h b/source/module_base/array_pool.h index 0d9d175c1d..5fde2c4176 100644 --- a/source/module_base/array_pool.h +++ b/source/module_base/array_pool.h @@ -40,7 +40,7 @@ namespace ModuleBase : nr(nr_in), nc(nc_in) { - this->ptr_1D = new T[nr * nc]; + this->ptr_1D = new T[nr * nc](); this->ptr_2D = new T*[nr]; for (int ir = 0; ir < nr; ++ir) this->ptr_2D[ir] = &this->ptr_1D[ir * nc]; diff --git a/source/module_base/test/vector3_test.cpp b/source/module_base/test/vector3_test.cpp index 0ac51c72f8..9a1cf49254 100644 --- a/source/module_base/test/vector3_test.cpp +++ b/source/module_base/test/vector3_test.cpp @@ -64,6 +64,9 @@ * - VneV * - overload operator "!=" to assert * - the inequality between two 3d vectors + * - VltV + * - overload operator "<" to sort + * - the "less than" relationship between two 3d vectors * - StdOutV * - overload operator "<<" to print out * - a 3d vectors on standard output @@ -703,6 +706,22 @@ TEST_F(Vector3Test,VneV) EXPECT_TRUE(wp != w); } +TEST_F(Vector3Test, VltV) +{ + ModuleBase::Vector3 u, up; + u.set(da, db, dc); + up.set(dc, db, da); + EXPECT_TRUE(u < up); + ModuleBase::Vector3 v, vp; + v.set(fa, fb, fc); + vp.set(fa, fb, fc); + EXPECT_FALSE(v < vp); + ModuleBase::Vector3 w, wp; + w.set(ia, ib, ic); + wp.set(ib, ib, ic); + EXPECT_TRUE(w < wp); +} + TEST_F(Vector3Test,StdOutV) { // double Vector3 diff --git a/source/module_base/vector3.h b/source/module_base/vector3.h index fda97c5482..3540018ae0 100644 --- a/source/module_base/vector3.h +++ b/source/module_base/vector3.h @@ -365,6 +365,22 @@ template inline Vector3 cross(const Vector3 &u, const Vector3 // (u.z * (v.x * w.y - v.y * w.x))); // } +// Overload the < operator for sorting +template bool operator<(const Vector3 &u, const Vector3 &v) +{ + if (u.x < v.x) + return true; + if (u.x > v.x) + return false; + if (u.y < v.y) + return true; + if (u.y > v.y) + return false; + if (u.z < v.z) + return true; + return false; +} + // whether m1 != m2 template inline bool operator!=(const Vector3 &u, const Vector3 &v) { diff --git a/source/module_basis/module_ao/ORB_atomic.h b/source/module_basis/module_ao/ORB_atomic.h index 66cd5fa8a6..2539e0a0c9 100644 --- a/source/module_basis/module_ao/ORB_atomic.h +++ b/source/module_basis/module_ao/ORB_atomic.h @@ -66,6 +66,7 @@ class Numerical_Orbital const inline Numerical_Orbital_Lm& PhiLN( const int &L, const int &N)const { + assert(this->phiLN != nullptr); return this->phiLN[ this->find_chi(L, N) ]; } diff --git a/source/module_elecstate/elecstate_lcao.cpp b/source/module_elecstate/elecstate_lcao.cpp index 748ef7a9b8..bae90aba79 100644 --- a/source/module_elecstate/elecstate_lcao.cpp +++ b/source/module_elecstate/elecstate_lcao.cpp @@ -8,7 +8,8 @@ #include "module_hamilt_lcao/module_gint/grid_technique.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_parameter/parameter.h" -#include "elecstate_lcao_cal_tau.h" + +#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h" #include @@ -59,13 +60,17 @@ void ElecStateLCAO>::psiToRho(const psi::Psigint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); this->gint_k->cal_gint(&inout); +#else + ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho); +#endif if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - elecstate::lcao_cal_tau_k(gint_k, this->charge); + this->cal_tau(psi); } this->charge->renormalize_rho(); @@ -92,15 +97,17 @@ void ElecStateLCAO::psiToRho(const psi::Psi& psi) //------------------------------------------------------------ ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!"); +#ifndef __NEW_GINT this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint - Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); - this->gint_gamma->cal_gint(&inout); +#else + ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho); +#endif if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { - elecstate::lcao_cal_tau_gamma(gint_gamma, this->charge); + this->cal_tau(psi); } this->charge->renormalize_rho(); @@ -158,17 +165,25 @@ void ElecStateLCAO::dmToRho(std::vector pexsi_DM, std::vectorgint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho, PARAM.inp.nspin); this->gint_gamma->cal_gint(&inout); +#else + ModuleGint::cal_gint_rho(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->rho); +#endif if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { for (int is = 0; is < PARAM.inp.nspin; is++) { ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx); } +#ifndef __NEW_GINT Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau); this->gint_gamma->cal_gint(&inout1); +#else + ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r); +#endif } this->charge->renormalize_rho(); diff --git a/source/module_elecstate/elecstate_lcao.h b/source/module_elecstate/elecstate_lcao.h index c85f6e27e5..4beeb017f0 100644 --- a/source/module_elecstate/elecstate_lcao.h +++ b/source/module_elecstate/elecstate_lcao.h @@ -46,6 +46,7 @@ class ElecStateLCAO : public ElecState // virtual void psiToRho(const psi::Psi& psi) override; // return current electronic density rho, as a input for constructing Hamiltonian // const double* getRho(int spin) const override; + virtual void cal_tau(const psi::Psi& psi) override; // update charge density for next scf step // void getNewRho() override; diff --git a/source/module_elecstate/elecstate_lcao_cal_tau.cpp b/source/module_elecstate/elecstate_lcao_cal_tau.cpp index c7d83bd1e9..9f2c0bd40b 100644 --- a/source/module_elecstate/elecstate_lcao_cal_tau.cpp +++ b/source/module_elecstate/elecstate_lcao_cal_tau.cpp @@ -1,56 +1,49 @@ #include "elecstate_lcao.h" -#include "elecstate_lcao_cal_tau.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h" + #include "module_base/timer.h" namespace elecstate { // calculate the kinetic energy density tau, multi-k case -void lcao_cal_tau_k(Gint_k* gint_k, - Charge* charge) +template <> +void ElecStateLCAO>::cal_tau(const psi::Psi>& psi) { ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); for (int is = 0; is < PARAM.inp.nspin; is++) { - ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx); + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); } - Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); - gint_k->cal_gint(&inout1); - +#ifndef __NEW_GINT + Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); + this->gint_k->cal_gint(&inout1); +#else + ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r); +#endif ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); return; } // calculate the kinetic energy density tau, gamma-only case -void lcao_cal_tau_gamma(Gint_Gamma* gint_gamma, - Charge* charge) +template <> +void ElecStateLCAO::cal_tau(const psi::Psi& psi) { ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); for (int is = 0; is < PARAM.inp.nspin; is++) { - ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx); + ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx); } - Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); - gint_gamma->cal_gint(&inout1); +#ifndef __NEW_GINT + Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau, PARAM.inp.nspin); + this->gint_gamma->cal_gint(&inout1); +#else + ModuleGint::cal_gint_tau(this->DM->get_DMR_vector(), PARAM.inp.nspin, this->charge->kin_r); +#endif ModuleBase::timer::tick("ElecStateLCAO", "cal_tau"); return; } -template <> -void lcao_cal_tau(Gint_Gamma* gint_gamma, - Gint_k* gint_k, - Charge* charge) -{ - lcao_cal_tau_gamma(gint_gamma, charge); -} -template <> -void lcao_cal_tau>(Gint_Gamma* gint_gamma, - Gint_k* gint_k, - Charge* charge) -{ - lcao_cal_tau_k(gint_k, charge); -} - -} // namespace elecstate \ No newline at end of file +} \ No newline at end of file diff --git a/source/module_elecstate/elecstate_lcao_cal_tau.h b/source/module_elecstate/elecstate_lcao_cal_tau.h deleted file mode 100644 index c0cfbc078a..0000000000 --- a/source/module_elecstate/elecstate_lcao_cal_tau.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef ELECSTATE_LCAO_CAL_TAU_H -#define ELECSTATE_LCAO_CAL_TAU_H -#include "module_elecstate/module_charge/charge.h" -#include "module_hamilt_lcao/module_gint/gint_gamma.h" -#include "module_hamilt_lcao/module_gint/gint_k.h" -namespace elecstate -{ - - void lcao_cal_tau_k(Gint_k* gint_k, - Charge* charge); - - void lcao_cal_tau_gamma(Gint_Gamma* gint_gamma, - Charge* charge); - - template - void lcao_cal_tau(Gint_Gamma* gint_gamma, - Gint_k* gint_k, - Charge* charge); - -} -#endif \ No newline at end of file diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 922d3e19bf..898cd55592 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -34,7 +34,6 @@ #include "module_base/global_function.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_elecstate/cal_ux.h" -#include "module_elecstate/elecstate_lcao_cal_tau.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need DeePKS_init @@ -933,7 +932,8 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) // 1) calculate the kinetic energy density tau, sunliang 2024-09-18 if (PARAM.inp.out_elf[0] > 0) { - elecstate::lcao_cal_tau(&(this->GG), &(this->GK), this->pelec->charge); + assert(this->psi != nullptr); + this->pelec->cal_tau(*(this->psi)); } //! 2) call after_scf() of ESolver_KS diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index ed2c180bed..8a620291f7 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -5,6 +5,8 @@ // for grid integration #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_info.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint.h" #ifdef __DEEPKS #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #endif diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 0bf61d6947..514b9f8981 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -115,6 +115,26 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) dpsi_u, d2psi_u, PARAM.inp.nstream); + +#ifdef __NEW_GINT + auto gint_info = std::make_shared( + this->pw_big->nbx, + this->pw_big->nby, + this->pw_big->nbz, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + 0, + 0, + this->pw_big->nbzp_start, + this->pw_big->nbx, + this->pw_big->nby, + this->pw_big->nbzp, + orb_.Phi, + ucell, + this->gd); + ModuleGint::Gint::init_gint_info(gint_info); +#endif psi_u.clear(); psi_u.shrink_to_fit(); dpsi_u.clear(); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp index aeb5d55c01..6e70b8c8fb 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp @@ -4,6 +4,7 @@ #include "module_base/tool_title.h" #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_cell/unitcell.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h" namespace hamilt { @@ -55,10 +56,53 @@ void Veff>::initialize_HR(const UnitCell* ucell_in, const G ModuleBase::timer::tick("Veff", "initialize_HR"); } +template<> +void Veff>::contributeHR() +{ + ModuleBase::TITLE("Veff", "contributeHR"); + ModuleBase::timer::tick("Veff", "contributeHR"); + //----------------------------------------- + //(1) prepare data for this k point. + // copy the local potential from array. + //----------------------------------------- + double* vr_eff1 = this->pot->get_effective_v(this->current_spin); + double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin); +#ifndef __NEW_GINT + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + { + Gint_inout inout(vr_eff1, vofk_eff1, Gint_Tools::job_type::vlocal_meta); + this->GG->cal_vlocal(&inout, this->new_e_iteration); + } + else + { + Gint_inout inout(vr_eff1, Gint_Tools::job_type::vlocal); + this->GG->cal_vlocal(&inout, this->new_e_iteration); + } + this->GG->transfer_pvpR(this->hR,this->ucell); + this->new_e_iteration = false; +#else + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + { + ModuleGint::cal_gint_vl_metagga(vr_eff1, vofk_eff1, this->hR); + } + else + { + ModuleGint::cal_gint_vl(vr_eff1, this->hR); + } +#endif -template -void Veff>::contributeHR() + if(this->nspin == 2) + { + this->current_spin = 1 - this->current_spin; + } + + ModuleBase::timer::tick("Veff", "contributeHR"); + return; +} + +template<> +void Veff, double>>::contributeHR() { ModuleBase::TITLE("Veff", "contributeHR"); ModuleBase::timer::tick("Veff", "contributeHR"); @@ -69,12 +113,7 @@ void Veff>::contributeHR() double* vr_eff1 = this->pot->get_effective_v(this->current_spin); double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin); - //-------------------------------------------- - //(2) check if we need to calculate - // pvpR = < phi0 | v(spin) | phiR> for a new spin. - //-------------------------------------------- - // GlobalV::ofs_running << " (spin change)" << std::endl; - +#ifndef __NEW_GINT // if you change the place of the following code, // rememeber to delete the #include if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) @@ -89,75 +128,82 @@ void Veff>::contributeHR() this->GK->cal_gint(&inout); } - // added by zhengdy-soc, for non-collinear case - // integral 4 times, is there any method to simplify? - if (this->nspin == 4) + this->GK->transfer_pvpR(this->hR,this->ucell,this->gd); +#else + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) { - for (int is = 1; is < 4; is++) - { - vr_eff1 = this->pot->get_effective_v(is); - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) - { - vofk_eff1 = this->pot->get_effective_vofk(is); - } - - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) - { - Gint_inout inout(vr_eff1, vofk_eff1, is, Gint_Tools::job_type::vlocal_meta); - this->GK->cal_gint(&inout); - } - else - { - Gint_inout inout(vr_eff1, is, Gint_Tools::job_type::vlocal); - this->GK->cal_gint(&inout); - } - } + ModuleGint::cal_gint_vl_metagga(vr_eff1, vofk_eff1, this->hR); } - this->GK->transfer_pvpR(this->hR,this->ucell,this->gd); + else + { + ModuleGint::cal_gint_vl(vr_eff1, this->hR); + } +#endif - if(this->nspin == 2) { this->current_spin = 1 - this->current_spin; -} + if(this->nspin == 2) + { + this->current_spin = 1 - this->current_spin; + } ModuleBase::timer::tick("Veff", "contributeHR"); return; } -// special case of gamma-only template<> -void Veff>::contributeHR(void) +void Veff, std::complex>>::contributeHR() { ModuleBase::TITLE("Veff", "contributeHR"); ModuleBase::timer::tick("Veff", "contributeHR"); - //----------------------------------------- - //(1) prepare data for this k point. - // copy the local potential from array. - //----------------------------------------- - const double* vr_eff1 = this->pot->get_effective_v(this->current_spin); - const double* vofk_eff1 = this->pot->get_effective_vofk(this->current_spin); - - //-------------------------------------------- - // (3) folding matrix, - // and diagonalize the H matrix (T+Vl+Vnl). - //-------------------------------------------- - - if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) +#ifndef __NEW_GINT + double* vr_eff1 = nullptr; + double* vofk_eff1 = nullptr; + for (int is = 0; is < 4; is++) { - Gint_inout inout(vr_eff1, vofk_eff1, Gint_Tools::job_type::vlocal_meta); - this->GG->cal_vlocal(&inout, this->new_e_iteration); + vr_eff1 = this->pot->get_effective_v(is); + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + { + vofk_eff1 = this->pot->get_effective_vofk(is); + } + + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + { + Gint_inout inout(vr_eff1, vofk_eff1, is, Gint_Tools::job_type::vlocal_meta); + this->GK->cal_gint(&inout); + } + else + { + Gint_inout inout(vr_eff1, is, Gint_Tools::job_type::vlocal); + this->GK->cal_gint(&inout); + } } - else + this->GK->transfer_pvpR(this->hR,this->ucell,this->gd); +#else + std::vector vr_eff(4, nullptr); + std::vector vofk_eff(4, nullptr); + for (int is = 0; is < 4; is++) { - Gint_inout inout(vr_eff1, Gint_Tools::job_type::vlocal); - this->GG->cal_vlocal(&inout, this->new_e_iteration); + vr_eff[is] = this->pot->get_effective_v(is); + if(XC_Functional::get_func_type()==3 || XC_Functional::get_func_type()==5) + { + vofk_eff[is] = this->pot->get_effective_vofk(is); + if(is == 3) + { + ModuleGint::cal_gint_vl_metagga(vr_eff, vofk_eff, this->hR); + } + } + else + { + if(is == 3) + { + ModuleGint::cal_gint_vl(vr_eff, this->hR); + } + } } - this->GG->transfer_pvpR(this->hR,this->ucell); - - this->new_e_iteration = false; - - if(this->nspin == 2) this->current_spin = 1 - this->current_spin; +#endif ModuleBase::timer::tick("Veff", "contributeHR"); + return; } // definition of class template should in the end of file to avoid compiling warning @@ -166,4 +212,4 @@ template class Veff>; template class Veff, double>>; template class Veff, std::complex>>; -} +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp index feb567eca2..e3cdcaf2ee 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress_gint.hpp @@ -3,6 +3,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h" #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_parameter/parameter.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h" namespace PulayForceStress { template @@ -17,8 +18,10 @@ 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; + +#ifndef __NEW_GINT + if (set_dmr_gint) { gint.transfer_DM2DtoGrid(dm.get_DMR_vector()); } // 2d block to grid for (int is = 0; is < nspin; ++is) { const double* vr_eff1 = pot->get_effective_v(is); @@ -35,6 +38,28 @@ namespace PulayForceStress gint.cal_gint(&inout); } } +#else + std::vector vr_eff(nspin, nullptr); + std::vector vofk_eff(nspin, nullptr); + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + for (int is = 0; is < nspin; ++is) + { + 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 + { + for(int is = 0; is < nspin; ++is) + { + vr_eff[is] = pot->get_effective_v(is); + } + ModuleGint::cal_gint_fvl(nspin, vr_eff, dm.get_DMR_vector(), isforce, isstress, &f, &s); + } +#endif + if (isstress) { StressTools::stress_fill(-1.0, ucell.omega, s); } } } \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/CMakeLists.txt b/source/module_hamilt_lcao/module_gint/CMakeLists.txt index e0e1e2e4cb..7b43114adb 100644 --- a/source/module_hamilt_lcao/module_gint/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_gint/CMakeLists.txt @@ -29,6 +29,31 @@ list(APPEND objects init_orb.cpp ) +if(NEW_GINT) + list(APPEND objects + temp_gint/biggrid_info.cpp + temp_gint/big_grid.cpp + temp_gint/divide_info.cpp + temp_gint/gint_atom.cpp + temp_gint/gint_info.cpp + temp_gint/gint.cpp + temp_gint/gint_vl.cpp + temp_gint/gint_vl_metagga.cpp + temp_gint/gint_vl_nspin4.cpp + temp_gint/gint_vl_metagga_nspin4.cpp + temp_gint/gint_rho.cpp + temp_gint/gint_tau.cpp + temp_gint/gint_fvl.cpp + temp_gint/gint_fvl_meta.cpp + temp_gint/localcell_info.cpp + temp_gint/phi_operator.cpp + temp_gint/set_ddphi.cpp + temp_gint/unitcell_info.cpp + temp_gint/gint_common.cpp + temp_gint/gint_interface.cpp + ) +endif() + if(USE_CUDA) list(APPEND objects gint_gpu_interface.cpp diff --git a/source/module_hamilt_lcao/module_gint/gint_force_cpu_interface.cpp b/source/module_hamilt_lcao/module_gint/gint_force_cpu_interface.cpp index 0222f0234b..0c08ac019f 100644 --- a/source/module_hamilt_lcao/module_gint/gint_force_cpu_interface.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_force_cpu_interface.cpp @@ -107,7 +107,7 @@ void Gint::gint_kernel_force(Gint_inout* inout) { //do integration to get stress this-> cal_meshball_stress(na_grid, block_index.data(), psir_vlbr3_DM.get_ptr_1D(), - dpsirr_ylm.get_ptr_1D(), svl_dphi_thread); + dpsirr_ylm.get_ptr_1D(), svl_dphi_thread); } } #pragma omp critical(gint) diff --git a/source/module_hamilt_lcao/module_gint/gint_tools.cpp b/source/module_hamilt_lcao/module_gint/gint_tools.cpp index 04bf875599..4979947505 100644 --- a/source/module_hamilt_lcao/module_gint/gint_tools.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_tools.cpp @@ -32,7 +32,8 @@ void get_vindex(const int bxyz, const int bx, const int by, const int bz, } } } - } + +} // here vindex refers to local potentials diff --git a/source/module_hamilt_lcao/module_gint/grid_bigcell.cpp b/source/module_hamilt_lcao/module_gint/grid_bigcell.cpp index 2db31ab94b..7f47218096 100644 --- a/source/module_hamilt_lcao/module_gint/grid_bigcell.cpp +++ b/source/module_hamilt_lcao/module_gint/grid_bigcell.cpp @@ -131,7 +131,6 @@ void Grid_BigCell::init_grid_expansion(const UnitCell& ucell,double* rcut) this->dxe = static_cast( this->orbital_rmax * g1) +1; this->dye = static_cast( this->orbital_rmax * g2) +1; this->dze = static_cast( this->orbital_rmax * g3) +1; - //xiaohui add 'PARAM.inp.out_level' line, 2015-09-16 if(PARAM.inp.out_level != "m") ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"extended fft grid",dxe,dye,dze); diff --git a/source/module_hamilt_lcao/module_gint/grid_technique.cpp b/source/module_hamilt_lcao/module_gint/grid_technique.cpp index 009b13e7ad..d97fde3a19 100644 --- a/source/module_hamilt_lcao/module_gint/grid_technique.cpp +++ b/source/module_hamilt_lcao/module_gint/grid_technique.cpp @@ -10,6 +10,8 @@ #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_hsolver/kernels/cuda/helper_cuda.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_helper.h" + Grid_Technique::Grid_Technique() { #if ((defined __CUDA) /* || (defined __ROCM) */) if (PARAM.inp.device == "gpu") { @@ -114,7 +116,7 @@ void Grid_Technique::set_pbc_grid(const int& ncx_in, this->init_meshball(); - this->init_atoms_on_grid(ny, nplane, startz_current, ucell); + this->init_atoms_on_grid(ny, nplane, ucell); this->init_ijr_and_nnrg(ucell, gd); this->cal_trace_lo(ucell); @@ -129,8 +131,7 @@ void Grid_Technique::set_pbc_grid(const int& ncx_in, } void Grid_Technique::get_startind(const int& ny, - const int& nplane, - const int& startz_current) { + const int& nplane) { ModuleBase::TITLE("Grid_Technique", "get_startind"); assert(nbxx >= 0); @@ -155,7 +156,7 @@ void Grid_Technique::get_startind(const int& ny, ix = ibx * this->bx; iy = iby * this->by; - iz = (ibz + nbzp_start) * this->bz - startz_current; + iz = ibz * this->bz; int ind = iz + iy * nplane + ix * ny * nplane; @@ -170,12 +171,11 @@ void Grid_Technique::get_startind(const int& ny, // mohan add 2021-04-06 void Grid_Technique::init_atoms_on_grid(const int& ny, const int& nplane, - const int& startz_current, const UnitCell& ucell) { ModuleBase::TITLE("Grid_Technique", "init_atoms_on_grid"); assert(nbxx >= 0); - this->get_startind(ny, nplane, startz_current); + this->get_startind(ny, nplane); // (1) prepare data. // counting the number of atoms whose orbitals have @@ -378,6 +378,8 @@ void Grid_Technique::init_atoms_on_grid2(const int* index2normal, int count = 0; this->how_many_atoms = std::vector(nbxx, 0); ModuleBase::Memory::record("GT::how many atoms", sizeof(int) * nbxx); + std::vector coord_x(total_atoms_on_grid* bxyz, 0.0); + std::vector coords3(bxyz * 3, 0.0); for(int iat = 0; iat < ucell.nat; iat++) { const int it = ucell.iat2it[iat]; @@ -431,12 +433,24 @@ void Grid_Technique::init_atoms_on_grid2(const int* index2normal, this->which_atom[index] = iat; this->which_bigcell[index] = im; this->which_unitcell[index] = index2ucell[extgrid]; + for(int imcell = 0; imcell < this -> bxyz; imcell++) + { + const double dr_x = this->meshcell_pos[imcell][0] + dr_x_part; + coord_x[index * bxyz + imcell] = dr_x; + } ++count; ++how_many_atoms[bcell_idx_on_proc]; } } } + for(int i = 0; i < this->bxyz; i++) + { + for(int j = 0; j < 3; j++) + { + coords3[i * 3 + j] = this->meshcell_pos[i][j]; + } + } assert(count == total_atoms_on_grid); return; } @@ -543,7 +557,6 @@ void Grid_Technique::init_ijr_and_nnrg(const UnitCell& ucell, const Grid_Driver& ModuleBase::TITLE("Grid_Technique", "init_ijr_and_nnrg"); hamilt::HContainer hRGint_tmp(ucell.nat); - // prepare the row_index and col_index for construct AtomPairs, they are // same, name as orb_index std::vector orb_index(ucell.nat + 1); diff --git a/source/module_hamilt_lcao/module_gint/grid_technique.h b/source/module_hamilt_lcao/module_gint/grid_technique.h index 6c3c26aba9..caa05cb310 100644 --- a/source/module_hamilt_lcao/module_gint/grid_technique.h +++ b/source/module_hamilt_lcao/module_gint/grid_technique.h @@ -133,7 +133,6 @@ class Grid_Technique : public Grid_MeshBall { // atoms on meshball void init_atoms_on_grid(const int& ny, const int& nplane, - const int& startz_current, const UnitCell& ucell); void init_atoms_on_grid2(const int* index2normal, const UnitCell& ucell); // initialize the ijr_info and nnrg @@ -142,8 +141,7 @@ class Grid_Technique : public Grid_MeshBall { void cal_trace_lo(const UnitCell& ucell); void check_bigcell(int* ind_bigcell, char* bigcell_on_processor); void get_startind(const int& ny, - const int& nplane, - const int& startz_current); + const int& nplane); #if ((defined __CUDA) /* || (defined __ROCM) */) public: diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.cpp new file mode 100644 index 0000000000..d972cd90bb --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.cpp @@ -0,0 +1,112 @@ +#include "big_grid.h" + +namespace ModuleGint +{ +std::shared_ptr BigGrid::localcell_info_ = nullptr; +std::shared_ptr BigGrid::unitcell_info_ = nullptr; +std::shared_ptr BigGrid::biggrid_info_ = nullptr; + +BigGrid::BigGrid(int idx): idx_(idx){} + +void BigGrid::add_atom(const GintAtom* atom) +{ + atoms_.push_back(atom); +} + +int BigGrid::get_mgrid_phi_len() const +{ + int len = 0; + for(const auto& atom : atoms_) + { + len += atom->get_nw(); + } + return len; +} + +void BigGrid::set_atoms_startidx(std::vector& startidx) const +{ + startidx.resize(atoms_.size()); + startidx[0] = 0; + for(int i = 1; i < atoms_.size(); ++i) + { + startidx[i] = startidx[i-1] + atoms_[i-1]->get_nw(); + } +} + +void BigGrid::set_atoms_phi_len(std::vector& phi_len) const +{ + phi_len.resize(atoms_.size()); + for(int i = 0; i < atoms_.size(); ++i) + { + phi_len[i] = atoms_[i]->get_nw(); + } +} + +void BigGrid::set_mgrids_coord(std::vector& coord) const +{ + coord.resize(biggrid_info_->get_mgrids_num()); + Vec3d this_bgrid_coord = localcell_info_->get_bgrid_global_coord_3D(idx_); + for(int im = 0; im < biggrid_info_->get_mgrids_num(); ++im) + { + coord[im] = biggrid_info_->get_mgrid_coord(im) + this_bgrid_coord; + } +} + +void BigGrid::set_mgrids_local_idx(std::vector& mgrids_idx) const +{ + auto index_3d = localcell_info_->bgrid_idx_1Dto3D(idx_); + Vec3i startidx( + index_3d.x * biggrid_info_->get_nmx(), + index_3d.y * biggrid_info_->get_nmy(), + index_3d.z * biggrid_info_->get_nmz()); + mgrids_idx.resize(0); + for(int ix = 0; ix < biggrid_info_->get_nmx(); ++ix) + { + for(int iy = 0; iy < biggrid_info_->get_nmy(); ++iy) + { + for(int iz = 0; iz < biggrid_info_->get_nmz(); ++iz) + { + Vec3i idx_3d(startidx.x + ix, startidx.y + iy, startidx.z + iz); + mgrids_idx.push_back(localcell_info_->mgrid_idx_3Dto1D(idx_3d)); + } + } + } +} + +void BigGrid::set_atom_relative_coords(const Vec3i bgrid_idx, const Vec3d tau_in_bgrid, std::vector& atom_coord) const +{ + Vec3i this_bgrid_idx = localcell_info_->get_bgrid_global_idx_3D(idx_); + + // the relative coordinates of this big grid and the atom + Vec3d bgrid_relative_coord + = unitcell_info_->get_relative_coord(bgrid_idx, this_bgrid_idx) + tau_in_bgrid; + + atom_coord.resize(biggrid_info_->get_mgrids_num()); + for(int im = 0; im < biggrid_info_->get_mgrids_num(); ++im) + { + const Vec3d& mcell_coord = biggrid_info_->get_mgrid_coord(im); + atom_coord[im] = mcell_coord - bgrid_relative_coord; + } +} + + +void BigGrid::set_atom_relative_coords(const GintAtom* atom, std::vector& atom_coord) const +{ + set_atom_relative_coords(atom->get_bgrid_idx(), atom->get_tau_in_bgrid(), atom_coord); +} + +bool BigGrid::is_atom_on_bgrid(const GintAtom* atom) const +{ + std::vector coords; + this->set_atom_relative_coords(atom, coords); + for(const auto& dist : coords) + { + if(dist.norm() <= atom->get_rcut()) + { + return true; + } + } + return false; +} + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.h b/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.h new file mode 100644 index 0000000000..c1d5596e13 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/big_grid.h @@ -0,0 +1,92 @@ +#pragma once + +#include +#include +#include "gint_type.h" +#include "biggrid_info.h" +#include "localcell_info.h" +#include "unitcell_info.h" +#include "gint_atom.h" + +namespace ModuleGint +{ + +class BigGrid +{ + public: + // constructor + BigGrid(int idx); + + static void init_localcell_info(std::shared_ptr localcell_info) { localcell_info_ = localcell_info; }; + static void init_unitcell_info(std::shared_ptr unitcell_info) { unitcell_info_ = unitcell_info; }; + static void init_bgrid_info(std::shared_ptr biggrid_info) { biggrid_info_ = biggrid_info; }; + + // getter functions + int get_idx() const { return idx_; }; + std::shared_ptr get_localcell_info() const { return localcell_info_; }; + std::shared_ptr get_unitcell_info() const {return unitcell_info_; }; + std::shared_ptr get_bgrid_info() const { return biggrid_info_; }; + const std::vector& get_atoms() const { return atoms_; }; + const GintAtom* get_atom(int i) const { return atoms_[i]; }; + + // get the number of meshgrids in the big grid + int get_mgrids_num() const { return biggrid_info_->get_mgrids_num(); }; + + // get the number of atoms that can affect the big grid + int get_atoms_num() const { return atoms_.size(); }; + + // add an atom to the big grid + void add_atom(const GintAtom* atom); + + // get the total number of phi of a meshgrid + // return: (\sum_{i=0}^{atoms_->size()} atoms_[i]->nw) + int get_mgrid_phi_len() const; + + // set the start index of the phi of each atom + // return: vector[i] = \sum_{j=0}^{i-1} atoms_[j]->nw + void set_atoms_startidx(std::vector& startidx) const; + + // set the length of phi of each atom + void set_atoms_phi_len(std::vector& phi_len) const; + + // set the coordinates of the meshgrids of the big grid + void set_mgrids_coord(std::vector& coord) const; + + // set the 1D index of the meshgrids in the local cell + void set_mgrids_local_idx(std::vector& mgrids_idx) const; + + /** + * @brief Set the coordinates of the meshgrids of the big grid relative to an atom + * + * @param bgrid_idx the 3D index of the big grid, which contains the atom, in the unitcell + * @param tau_in_bgrid the cartesian coordinate of the atom relative to the big grid containing it + * @param atom_coord the relative cartesian coordinates of the atom and the meshgrids + */ + void set_atom_relative_coords(const Vec3i bgrid_idx, const Vec3d tau_in_bgrid, std::vector& atom_coord) const; + + // a wrapper function to get the relative coordinates of the atom and the meshgrids + void set_atom_relative_coords(const GintAtom* atom, std::vector& atom_coord) const; + + // if the atom affects the big grid, return true, otherwise false + // note when we say an atom affects a big grid, it does not mean that the atom affects all the meshgrid on the big grid, + // it may only affect a part of them. + bool is_atom_on_bgrid(const GintAtom* atom) const; + + private: + // atoms that can affect the big grid + std::vector atoms_; + + // the 1D index of the big grid in the local cell + const int idx_; + + // local cell info + static std::shared_ptr localcell_info_; + + // unitcell info + static std::shared_ptr unitcell_info_; + + // the big grid info + static std::shared_ptr biggrid_info_; +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.cpp new file mode 100644 index 0000000000..f69a7ef0ec --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.cpp @@ -0,0 +1,66 @@ +#include "biggrid_info.h" +#include "gint_helper.h" +#include "gint_type.h" + +namespace ModuleGint +{ + +BigGridInfo::BigGridInfo( + Vec3d biggrid_vec1, + Vec3d biggrid_vec2, + Vec3d biggrid_vec3, + int nmx, int nmy, int nmz) + : biggrid_vec1_(biggrid_vec1), + biggrid_vec2_(biggrid_vec2), + biggrid_vec3_(biggrid_vec3), + nmx_(nmx), nmy_(nmy), nmz_(nmz), nmxyz_(nmx*nmy*nmz) + { + // initialize the biggrid_latvec0_ + biggrid_latvec0_.e11 = biggrid_vec1_.x; + biggrid_latvec0_.e12 = biggrid_vec1_.y; + biggrid_latvec0_.e13 = biggrid_vec1_.z; + + biggrid_latvec0_.e21 = biggrid_vec2_.x; + biggrid_latvec0_.e22 = biggrid_vec2_.y; + biggrid_latvec0_.e23 = biggrid_vec2_.z; + + biggrid_latvec0_.e31 = biggrid_vec3_.x; + biggrid_latvec0_.e32 = biggrid_vec3_.y; + biggrid_latvec0_.e33 = biggrid_vec3_.z; + + // initialize the GT matrix + biggrid_GT_ = biggrid_latvec0_.Inverse(); + + // initialize the meshgrid_info_ + meshgrid_info_ = std::make_shared( + biggrid_vec1_ / static_cast(nmx), + biggrid_vec2_ / static_cast(nmy), + biggrid_vec3_ / static_cast(nmz)); + + // initialize the meshgrid_coords_ + meshgrid_coords_.resize(nmxyz_); + for(int index_1d = 0; index_1d < nmxyz_; index_1d++) + { + meshgrid_coords_[index_1d] = + meshgrid_info_->get_cartesian_coord(mgrid_idx_1Dto3D(index_1d)); + } + } + + Vec3i BigGridInfo::max_ext_bgrid_num(double r) const + { + const double g1 = sqrt(biggrid_GT_.e11 * biggrid_GT_.e11 + + biggrid_GT_.e21 * biggrid_GT_.e21 + + biggrid_GT_.e31 * biggrid_GT_.e31); + const double g2 = sqrt(biggrid_GT_.e12 * biggrid_GT_.e12 + + biggrid_GT_.e22 * biggrid_GT_.e22 + + biggrid_GT_.e32 * biggrid_GT_.e32); + const double g3 = sqrt(biggrid_GT_.e13 * biggrid_GT_.e13 + + biggrid_GT_.e23 * biggrid_GT_.e23 + + biggrid_GT_.e33 * biggrid_GT_.e33); + int ext_x = static_cast(r * g1) + 1; + int ext_y = static_cast(r * g2) + 1; + int ext_z = static_cast(r * g3) + 1; + return Vec3i(ext_x, ext_y, ext_z); + } + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.h new file mode 100644 index 0000000000..f8bcb79665 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/biggrid_info.h @@ -0,0 +1,99 @@ +#pragma once + +#include +#include "gint_type.h" +#include "gint_helper.h" +#include "meshgrid_info.h" + +namespace ModuleGint +{ + +/** + * @class BigGridInfo + * @brief This class stores some basic properties common to all big grids. + */ +class BigGridInfo +{ + public: + // constructor + BigGridInfo( + Vec3d biggrid_vec1, + Vec3d biggrid_vec2, + Vec3d biggrid_vec3, + int nmx, int nmy, int nmz); + + Vec3d get_cartesian_coord(const Vec3d& index_3d) const { return index_3d * biggrid_latvec0_; }; + Vec3d get_cartesian_coord(const Vec3i& index_3d) const { return index_3d * biggrid_latvec0_; }; + const Vec3d get_direct_coord(const Vec3d& cart_coord) const { return cart_coord * biggrid_GT_; }; + + // Return the maximum number of big grids that can fit inside a sphere of radius r, + // along the three lattice vector directions. + Vec3i max_ext_bgrid_num(double r) const; + + // get number of meshgrids along three lattice directions + int get_nmx() const { return nmx_; }; + int get_nmy() const { return nmy_; }; + int get_nmz() const { return nmz_; }; + int get_mgrids_num() const { return nmxyz_; }; + + const std::vector& get_mgrids_coord() const { return meshgrid_coords_; }; + const Vec3d& get_mgrid_coord(int index_1d) const { return meshgrid_coords_[index_1d]; }; + + std::shared_ptr get_mgrid_info() const { return meshgrid_info_; }; + + // get the 3D index of a meshgrid in the big grid from the 1D index + Vec3i mgrid_idx_1Dto3D(int index_1d) const + { + return index1Dto3D(index_1d, nmx_, nmy_, nmz_); + }; + + // get the 1D index of a meshgrid in the big grid from the 3D index + int mgrid_idx_3Dto1D(const Vec3i index_3d) const + { + return index3Dto1D(index_3d.x, index_3d.y, index_3d.z, nmx_, nmy_, nmz_); + }; + + private: + // basis vectors of the big grid + Vec3d biggrid_vec1_; + Vec3d biggrid_vec2_; + Vec3d biggrid_vec3_; + + // used to convert the (i, j, k) index of the big grid to the Cartesian coordinate + // if biggrid_vec1_ is row vector, + // then biggrid_latvec0_ = [biggrid_vec1_; biggrid_vec2_; biggrid_vec3_], + // (i, j, k) * biggrid_latvec0_ = (x, y, z) + Matrix3 biggrid_latvec0_; + + // used to convert the Cartesian coordinate to the (i, j, k) index of the big grid + // biggrid_GT_ = biggrid_latvec0_.Inverse() + // (x, y, z) * biggrid_GT_ = (i, j, k) + Matrix3 biggrid_GT_; + + //====================================================== + // some member variables related to meshgrid + //====================================================== + + // basic attributes of meshgrid + std::shared_ptr meshgrid_info_; + + // the number of meshgrids of a biggrid along the first basis vector + // nmx may be a confusing name, because it is not the number of meshgrids along x axis + // but it's used in the original code, so I keep it, maybe it will be changed later + int nmx_; + + // the number of meshgrids of a biggrid along the second basis vector + int nmy_; + + // the number of meshgrids of a biggrid along the third basis vector + int nmz_; + + // total number of meshgrids in the biggrid + int nmxyz_; + + // store the relative Cartesian coordinates of all meshgrids in the biggrid + // the size of vector is nbxyz_ + std::vector meshgrid_coords_; +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.cpp new file mode 100644 index 0000000000..c4a5b2a738 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.cpp @@ -0,0 +1,24 @@ +#include "divide_info.h" + +namespace ModuleGint +{ + +DivideInfo::DivideInfo( + int startidx_bx_old, int startidx_by_old, int startidx_bz_old, + int nbx_old, int nby_old, int nbz_old, + std::shared_ptr unitcell_info, bool is_redevided) + : startidx_bx_old_(startidx_bx_old), startidx_by_old_(startidx_by_old), startidx_bz_old_(startidx_bz_old), + nbx_old_(nbx_old), nby_old_(nby_old), nbz_old_(nbz_old), + startidx_bx_new_(startidx_bx_old), startidx_by_new_(startidx_by_old), startidx_bz_new_(startidx_bz_old), + nbx_new_(nbx_old), nby_new_(nby_old), nbz_new_(nbz_old), + unitcell_info_(unitcell_info), is_redivided_(is_redevided) + { + if(!is_redivided_) + { + localcell_info_ = std::make_shared(startidx_bx_new_, startidx_by_new_, startidx_bz_new_, + nbx_new_, nby_new_, nbz_new_, unitcell_info_); + } + // TODO: "implement the redivide function"; + } + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.h new file mode 100644 index 0000000000..374831eeef --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/divide_info.h @@ -0,0 +1,54 @@ +#pragma once + +#include "biggrid_info.h" +#include "unitcell_info.h" +#include "localcell_info.h" + +namespace ModuleGint +{ + +class DivideInfo +{ + public: + // constructor + DivideInfo( + int startidx_bx_old, int startidx_by_old, int startidx_bz_old, + int nbx_old, int nby_old, int nbz_old, + std::shared_ptr unitcell_info, bool is_redivided = false); + + // getter functions + std::shared_ptr get_localcell_info() const { return localcell_info_; } + bool get_is_redivided() const { return is_redivided_; } + + private: + // if the grid is redivided, is_redeiided_ is true + bool is_redivided_; + + // the old start index of the local cell + int startidx_bx_old_; + int startidx_by_old_; + int startidx_bz_old_; + + // the old number of big grids in the local cell + int nbx_old_; + int nby_old_; + int nbz_old_; + + // the new start index of the local cell + int startidx_bx_new_; + int startidx_by_new_; + int startidx_bz_new_; + + // the new number of big grids in the local cell + int nbx_new_; + int nby_new_; + int nbz_new_; + + // the unitcell info + std::shared_ptr unitcell_info_; + + // the localcell info + std::shared_ptr localcell_info_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint.cpp new file mode 100644 index 0000000000..e766c46d9f --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint.cpp @@ -0,0 +1,8 @@ +#include "gint.h" + +namespace ModuleGint +{ + +std::shared_ptr Gint::gint_info_ = nullptr; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint.h new file mode 100644 index 0000000000..a14f014a6c --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint.h @@ -0,0 +1,28 @@ +#pragma once +#include +#include "gint_info.h" +#include "gint_type.h" + +namespace ModuleGint +{ + +class Gint +{ + public: + Gint() = default; + virtual ~Gint() = default; + + virtual void cal_gint() = 0; + + // note that gint_info_ is a static member variable + // it is shared by all instances of Gint + static void init_gint_info(std::shared_ptr gint_info) + { + gint_info_ = gint_info; + } + + protected: + static std::shared_ptr gint_info_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.cpp new file mode 100644 index 0000000000..2c8252092e --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.cpp @@ -0,0 +1,207 @@ +#include "module_base/ylm.h" +#include "module_base/array_pool.h" +#include "gint_atom.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +template +void GintAtom::set_phi(const std::vector& coords, const int stride, T* phi) const +{ + const int num_mgrids = coords.size(); + + // orb_ does not have the member variable dr_uniform + const double dr_uniform = orb_->PhiLN(0, 0).dr_uniform; + + // store the pointer to reduce repeated address fetching + std::vector p_psi_uniform(atom_->nw); + std::vector p_dpsi_uniform(atom_->nw); + for(int iw = 0; iw < atom_->nw; iw++) + { + if(atom_->iw2_new[iw]) + { + int l = atom_->iw2l[iw]; + int n = atom_->iw2n[iw]; + p_psi_uniform[iw] = orb_->PhiLN(l, n).psi_uniform.data(); + p_dpsi_uniform[iw] = orb_->PhiLN(l, n).dpsi_uniform.data(); + } + } + + // store the spherical harmonics + // it's outside the loop to reduce the vector allocation overhead + std::vector ylma; + + for(int im = 0; im < num_mgrids; im++) + { + const Vec3d& coord = coords[im]; + + // 1e-9 is to avoid division by zero + const double dist = coord.norm() < 1e-9 ? 1e-9 : coord.norm(); + if(dist > orb_->getRcut()) + { + // if the distance is larger than the cutoff radius, + // the wave function values are all zeros + ModuleBase::GlobalFunc::ZEROS(phi + im * stride, atom_->nw); + } + else + { + // spherical harmonics + // TODO: vectorize the sph_harm function, + // the vectorized function can be called once for all meshgrids in a biggrid + ModuleBase::Ylm::sph_harm(atom_->nwl, coord.x/dist, coord.y/dist, coord.z/dist, ylma); + // interpolation + + // these parameters are related to interpolation + // because once the distance from atom to grid point is known, + // we can obtain the parameters for interpolation and + // store them first! these operations can save lots of efforts. + const double position = dist / dr_uniform; + const int ip = static_cast(position); + const double dx = position - ip; + const double dx2 = dx * dx; + const double dx3 = dx2 * dx; + + const double c3 = 3.0 * dx2 - 2.0 * dx3; + const double c1 = 1.0 - c3; + const double c2 = (dx - 2.0 * dx2 + dx3) * dr_uniform; + const double c4 = (dx3 - dx2) * dr_uniform; + + // I'm not sure if the variable name 'psi' is appropriate + double psi = 0; + + for(int iw = 0; iw < atom_->nw; iw++) + { + if(atom_->iw2_new[iw]) + { + auto psi_uniform = p_psi_uniform[iw]; + auto dpsi_uniform = p_dpsi_uniform[iw]; + psi = c1 * psi_uniform[ip] + c2 * dpsi_uniform[ip] + + c3 * psi_uniform[ip + 1] + c4 * dpsi_uniform[ip + 1]; + } + phi[im * stride + iw] = psi * ylma[atom_->iw2_ylm[iw]]; + } + } + } +} + +template +void GintAtom::set_phi_dphi( + const std::vector& coords, const int stride, + T* phi, T* dphi_x, T* dphi_y, T* dphi_z) const +{ + const int num_mgrids = coords.size(); + + // orb_ does not have the member variable dr_uniform + const double dr_uniform = orb_->PhiLN(0, 0).dr_uniform; + + // store the pointer to reduce repeated address fetching + std::vector p_psi_uniform(atom_->nw); + std::vector p_dpsi_uniform(atom_->nw); + std::vector phi_nr_uniform(atom_->nw); + for (int iw=0; iw< atom_->nw; ++iw) + { + if ( atom_->iw2_new[iw] ) + { + int l = atom_->iw2l[iw]; + int n = atom_->iw2n[iw]; + p_psi_uniform[iw] = orb_->PhiLN(l, n).psi_uniform.data(); + p_dpsi_uniform[iw] = orb_->PhiLN(l, n).dpsi_uniform.data(); + phi_nr_uniform[iw] = orb_->PhiLN(l, n).nr_uniform; + } + } + + std::vector rly(std::pow(atom_->nwl + 1, 2)); + // TODO: replace array_pool with std::vector + ModuleBase::Array_Pool grly(std::pow(atom_->nwl + 1, 2), 3); + + for(int im = 0; im < num_mgrids; im++) + { + const Vec3d& coord = coords[im]; + // 1e-9 is to avoid division by zero + const double dist = coord.norm() < 1e-9 ? 1e-9 : coord.norm(); + + if(dist > orb_->getRcut()) + { + // if the distance is larger than the cutoff radius, + // the wave function values are all zeros + if(phi != nullptr) + { + ModuleBase::GlobalFunc::ZEROS(phi + im * stride, atom_->nw); + } + ModuleBase::GlobalFunc::ZEROS(dphi_x + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(dphi_y + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(dphi_z + im * stride, atom_->nw); + } + else + { + // spherical harmonics + // TODO: vectorize the sph_harm function, + // the vectorized function can be called once for all meshgrids in a biggrid + ModuleBase::Ylm::grad_rl_sph_harm(atom_->nwl, coord.x, coord.y, coord.z, rly.data(), grly.get_ptr_2D()); + + // interpolation + const double position = dist / dr_uniform; + const int ip = static_cast(position); + const double x0 = position - ip; + const double x1 = 1.0 - x0; + const double x2 = 2.0 - x0; + const double x3 = 3.0 - x0; + const double x12 = x1 * x2 / 6; + const double x03 = x0 * x3 / 2; + + double tmp, dtmp; + for(int iw = 0; iw < atom_->nw; ++iw) + { + // this is a new 'l', we need 1D orbital wave + // function from interpolation method. + if(atom_->iw2_new[iw]) + { + auto psi_uniform = p_psi_uniform[iw]; + auto dpsi_uniform = p_dpsi_uniform[iw]; + + if(ip >= phi_nr_uniform[iw] - 4) + { + tmp = dtmp = 0.0; + } + else + { + // use Polynomia Interpolation method to get the + // wave functions + + tmp = x12 * (psi_uniform[ip] * x3 + psi_uniform[ip + 3] * x0) + + x03 * (psi_uniform[ip + 1] * x2 - psi_uniform[ip + 2] * x1); + + dtmp = x12 * (dpsi_uniform[ip] * x3 + dpsi_uniform[ip + 3] * x0) + + x03 * (dpsi_uniform[ip + 1] * x2 - dpsi_uniform[ip + 2] * x1); + } + } // new l is used. + + // get the 'l' of this localized wave function + const int ll = atom_->iw2l[iw]; + const int idx_lm = atom_->iw2_ylm[iw]; + + const double rl = pow_int(dist, ll); + const double tmprl = tmp / rl; + + // 3D wave functions + if(phi != nullptr) + { + phi[im * stride + iw] = tmprl * rly[idx_lm]; + } + + // derivative of wave functions with respect to atom positions. + const double tmpdphi_rly = (dtmp - tmp * ll / dist) / rl * rly[idx_lm] / dist; + + dphi_x[im * stride + iw] = tmpdphi_rly * coord.x + tmprl * grly[idx_lm][0]; + dphi_y[im * stride + iw] = tmpdphi_rly * coord.y + tmprl * grly[idx_lm][1]; + dphi_z[im * stride + iw] = tmpdphi_rly * coord.z + tmprl * grly[idx_lm][2]; + } + } + } +} + +// explicit instantiation +template void GintAtom::set_phi(const std::vector& coords, const int stride, double* phi) const; +template void GintAtom::set_phi_dphi(const std::vector& coords, const int stride, double* phi, double* dphi_x, double* dphi_y, double* dphi_z) const; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h new file mode 100644 index 0000000000..660aa0ceab --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_atom.h @@ -0,0 +1,118 @@ +#pragma once + +#include "module_cell/atom_spec.h" +#include "module_basis/module_ao/ORB_atomic.h" +#include "gint_type.h" + +namespace ModuleGint +{ + +class GintAtom +{ + public: + // constructor + GintAtom( + const Atom* atom, + int ia, + int iat, + Vec3i biggrid_idx, + Vec3i unitcell_idx, + Vec3d tau_in_biggrid, + const Numerical_Orbital* orb) + : atom_(atom), ia_(ia), iat_(iat), biggrid_idx_(biggrid_idx), + unitcell_idx_(unitcell_idx), tau_in_biggrid_(tau_in_biggrid), + orb_(orb) {}; + + // getter functions + const Atom* get_atom() const { return atom_; }; + const int get_ia() const { return ia_; }; + const int get_iat() const { return iat_; }; + const Vec3i& get_bgrid_idx() const { return biggrid_idx_; }; + const Vec3i& get_unitcell_idx() const { return unitcell_idx_; }; + const Vec3i& get_R() const { return unitcell_idx_; }; + const Vec3d& get_tau_in_bgrid() const { return tau_in_biggrid_; }; + const Numerical_Orbital* get_orb() const { return orb_; }; + + const int get_nw() const { return atom_->nw; }; + const double get_rcut() const { return orb_->getRcut(); }; + + /** + * @brief Get the wave function values of the atom at a meshgrid. + * + * phi[(n-1)*stride] ~ phi[(n-1)*stride + nw] store the wave function values of the first atom at the nth meshgrid + * + * @param coords the cartesian coordinates of the meshgrids of a biggrid relative to the atom + * @param stride the stride of the phi array between two adjacent meshgrids + * @param phi array to store the wave function values + */ + template + void set_phi(const std::vector& coords, const int stride, T* phi) const; + + /** + * @brief Get the wave function values and its derivative + * + * The reason for combining the functions to solve the wave function values + * and wave function derivatives into one function is to improve efficiency. + * phi[(n-1)*stride] ~ phi[(n-1)*stride + nw] store the wave function values of the first atom at the nth meshgrid + * + * @param coords the cartesian coordinates of the meshgrids of a biggrid relative to the atom + * @param stride the stride of the phi array between two adjacent meshgrids + * @param phi array to store the wave function values + * @param dphi_x array to store the derivative wave functions in x direction + * @param dphi_y array to store the derivative wave functions in y direction + * @param dphi_z array to store the derivative wave functions in z direction + */ + template + void set_phi_dphi( + const std::vector& coords, const int stride, + T* phi, T* dphi_x, T* dphi_y, T* dphi_z) const; + + /** + * @brief Get the wave function values and its second derivative + * + * ddphi[(n-1)*stride] ~ ddphi[(n-1)*stride + nw] store the second derivative of + * wave function values of the atom at the first meshgrid + * + * @param coords the cartesian coordinates of the meshgrids of a biggrid relative to the atom + * @param stride the stride of the phi array between two adjacent meshgrids + * @param ddphi_xx array to store the second derivative wave functions in xx direction + * @param ddphi_xy array to store the second derivative wave functions in xy direction + * @param ddphi_xz array to store the second derivative wave functions in xz direction + * @param ddphi_yy array to store the second derivative wave functions in yy direction + * @param ddphi_yz array to store the second derivative wave functions in yz direction + * @param ddphi_zz array to store the second derivative wave functions in zz direction + */ + template + void set_ddphi( + const std::vector& coords, const int stride, + T* ddphi_xx, T* ddphi_xy, T* ddphi_xz, + T* ddphi_yy, T* ddphi_yz, T* ddphi_zz) const; + + private: + // the atom object + const Atom* atom_; + + // the global index of the atom + int iat_; + + // the global index of the atom among the same type of atoms + int ia_; + + // the index of big grid which contains this atom + Vec3i biggrid_idx_; + + // the index of the unitcell which contains this atom + Vec3i unitcell_idx_; + + // the relative Cartesian coordinates of this atom + // with respect to the big grid that contains it + Vec3d tau_in_biggrid_; + + // the numerical orbitals of this atom + // In fact, I think the Numerical_Orbital class + // should be a member of the Atom class, not the GintAtom class + const Numerical_Orbital* orb_; + +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.cpp new file mode 100644 index 0000000000..0bf8ad7726 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.cpp @@ -0,0 +1,196 @@ +#include "gint_common.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +#include "module_parameter/parameter.h" + +namespace ModuleGint +{ + +void compose_hRGint(std::shared_ptr> hRGint) +{ + for (int iap = 0; iap < hRGint->size_atom_pairs(); iap++) + { + auto& ap = hRGint->get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + if (iat1 > iat2) + { + // fill lower triangle matrix with upper triangle matrix + // the upper is + const hamilt::AtomPair* upper_ap = hRGint->find_pair(iat2, iat1); + const hamilt::AtomPair* lower_ap = hRGint->find_pair(iat1, iat2); +#ifdef __DEBUG + assert(upper_ap != nullptr); +#endif + for (int ir = 0; ir < ap.get_R_size(); ir++) + { + auto R_index = ap.get_R_index(ir); + auto upper_mat = upper_ap->find_matrix(-R_index); + auto lower_mat = lower_ap->find_matrix(R_index); + for (int irow = 0; irow < upper_mat->get_row_size(); ++irow) + { + for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) + { + lower_mat->get_value(icol, irow) = upper_ap->get_value(irow, icol); + } + } + } + } + } +} + +void compose_hRGint(std::vector>> hRGint_part, + std::shared_ptr>> hRGint_full) +{ + for (int iap = 0; iap < hRGint_full->size_atom_pairs(); iap++) + { + auto* ap = &hRGint_full->get_atom_pair(iap); + const int iat1 = ap->get_atom_i(); + const int iat2 = ap->get_atom_j(); + if (iat1 <= iat2) + { + hamilt::AtomPair>* upper_ap = ap; + hamilt::AtomPair>* lower_ap = hRGint_full->find_pair(iat2, iat1); + const hamilt::AtomPair* ap_nspin_0 = hRGint_part[0]->find_pair(iat1, iat2); + const hamilt::AtomPair* ap_nspin_3 = hRGint_part[3]->find_pair(iat1, iat2); + for (int ir = 0; ir < upper_ap->get_R_size(); ir++) + { + const auto R_index = upper_ap->get_R_index(ir); + auto upper_mat = upper_ap->find_matrix(R_index); + auto mat_nspin_0 = ap_nspin_0->find_matrix(R_index); + auto mat_nspin_3 = ap_nspin_3->find_matrix(R_index); + + // The row size and the col size of upper_matrix is double that of matrix_nspin_0 + for (int irow = 0; irow < mat_nspin_0->get_row_size(); ++irow) + { + for (int icol = 0; icol < mat_nspin_0->get_col_size(); ++icol) + { + upper_mat->get_value(2*irow, 2*icol) = mat_nspin_0->get_value(irow, icol) + mat_nspin_3->get_value(irow, icol); + upper_mat->get_value(2*irow+1, 2*icol+1) = mat_nspin_0->get_value(irow, icol) - mat_nspin_3->get_value(irow, icol); + } + } + + if (PARAM.globalv.domag) + { + const hamilt::AtomPair* ap_nspin_1 = hRGint_part[1]->find_pair(iat1, iat2); + const hamilt::AtomPair* ap_nspin_2 = hRGint_part[2]->find_pair(iat1, iat2); + const auto mat_nspin_1 = ap_nspin_1->find_matrix(R_index); + const auto mat_nspin_2 = ap_nspin_2->find_matrix(R_index); + for (int irow = 0; irow < mat_nspin_1->get_row_size(); ++irow) + { + for (int icol = 0; icol < mat_nspin_1->get_col_size(); ++icol) + { + upper_mat->get_value(2*irow, 2*icol+1) = mat_nspin_1->get_value(irow, icol) + std::complex(0.0, 1.0) * mat_nspin_2->get_value(irow, icol); + upper_mat->get_value(2*irow+1, 2*icol) = mat_nspin_1->get_value(irow, icol) - std::complex(0.0, 1.0) * mat_nspin_2->get_value(irow, icol); + } + } + } + + // fill the lower triangle matrix + if (iat1 < iat2) + { + auto lower_mat = lower_ap->find_matrix(-R_index); + for (int irow = 0; irow < upper_mat->get_row_size(); ++irow) + { + for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) + { + lower_mat->get_value(icol, irow) = conj(upper_mat->get_value(irow, icol)); + } + } + } + } + } + } +} + +template +void transfer_hRGint_to_hR(std::shared_ptr> hRGint, HContainer* hR) +{ +#ifdef __MPI + int size = 0; + MPI_Comm_size(MPI_COMM_WORLD, &size); + if (size == 1) + { + hR->add(*hRGint); + } + else + { + hamilt::transferSerials2Parallels(*hRGint, hR); + } +#else + hR->add(*hRGint); +#endif +} + +// gint_info should not have been a parameter, but it was added to initialize DMRGint_full +// In the future, we might try to remove the gint_info parameter +void transfer_DM_to_DMGint( + std::shared_ptr gint_info, + std::vector*> DM, + std::vector>> DMRGint) +{ + // To check whether input parameter DM2D has been initialized +#ifdef __DEBUG + assert(PARAM.inp.nspin == DM.size() + && "The size of DM should be equal to the number of spins!"); +#endif + + if (PARAM.inp.nspin != 4) + { + for (int is = 0; is < PARAM.inp.nspin; is++) + { +#ifdef __MPI + hamilt::transferParallels2Serials(*DM[is], DMRGint[is].get()); +#else + DMRGint[is]->set_zero(); + DMRGint[is]->add(*DM[is]); +#endif + } + } else // NSPIN=4 case + { +#ifdef __MPI + const int npol = 2; + std::shared_ptr> DM_full = gint_info->get_hr(npol); + hamilt::transferParallels2Serials(*DM[0], DM_full.get()); +#else + HContainer* DM_full = DM[0]; +#endif + std::vector tmp_pointer(4, nullptr); + for (int iap = 0; iap < DM_full->size_atom_pairs(); iap++) + { + auto& ap = DM_full->get_atom_pair(iap); + const int iat1 = ap.get_atom_i(); + const int iat2 = ap.get_atom_j(); + for (int ir = 0; ir < ap.get_R_size(); ir++) + { + const ModuleBase::Vector3 r_index = ap.get_R_index(ir); + for (int is = 0; is < 4; is++) + { + tmp_pointer[is] = + DMRGint[is]->find_matrix(iat1, iat2, r_index)->get_pointer(); + } + double* data_full = ap.get_pointer(ir); + for (int irow = 0; irow < ap.get_row_size(); irow += 2) + { + for (int icol = 0; icol < ap.get_col_size(); icol += 2) + { + *(tmp_pointer[0])++ = data_full[icol]; + *(tmp_pointer[1])++ = data_full[icol + 1]; + } + data_full += ap.get_col_size(); + for (int icol = 0; icol < ap.get_col_size(); icol += 2) + { + *(tmp_pointer[2])++ = data_full[icol]; + *(tmp_pointer[3])++ = data_full[icol + 1]; + } + data_full += ap.get_col_size(); + } + } + } + } +} + + +template void transfer_hRGint_to_hR(std::shared_ptr> hRGint, HContainer* hR); +template void transfer_hRGint_to_hR(std::shared_ptr>> hRGint, HContainer>* hR); +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.h new file mode 100644 index 0000000000..aa47fdf351 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_common.h @@ -0,0 +1,21 @@ +#pragma once +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_hamilt_lcao/module_gint/temp_gint/gint_info.h" + +namespace ModuleGint +{ + // fill the lower triangle matrix with the upper triangle matrix + void compose_hRGint(std::shared_ptr> hRGint); + // for nspin=4 case + void compose_hRGint(std::vector>> hRGint_part, + std::shared_ptr>> hRGint_full); + + template + void transfer_hRGint_to_hR(std::shared_ptr> hRGint, HContainer* hR); + + void transfer_DM_to_DMGint( + std::shared_ptr gint_info, + std::vector*> DM, + std::vector>> DMRGint); + +} diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.cpp new file mode 100644 index 0000000000..adbffddf91 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.cpp @@ -0,0 +1,95 @@ +#include "module_base/global_function.h" +#include "gint_fvl.h" +#include "gint_common.h" +#include "phi_operator.h" + +namespace ModuleGint +{ + +void Gint_fvl::cal_gint() +{ + init_DMRGint_(); + transfer_DM_to_DMGint(gint_info_, DMR_vec_, DMRGint_vec_); + cal_fvl_svl_(); +} + +void Gint_fvl::init_DMRGint_() +{ + DMRGint_vec_.resize(nspin_); + for (int is = 0; is < nspin_; is++) + { + DMRGint_vec_[is] = gint_info_->get_hr(); + } +} + +void Gint_fvl::cal_fvl_svl_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + std::vector phi_vldr3_DM; + std::vector dphi_x; + std::vector dphi_y; + std::vector dphi_z; + 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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + phi_vldr3_DM.resize(phi_len); + dphi_x.resize(phi_len); + dphi_y.resize(phi_len); + dphi_z.resize(phi_len); + phi_op.set_phi_dphi(phi.data(), dphi_x.data(), dphi_y.data(), dphi_z.data()); + for (int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_vldr3(vr_eff_[is], dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_dm(phi_vldr3.data(), *DMRGint_vec_[is], false, phi_vldr3_DM.data()); + if(isforce_) + { + phi_op.phi_dot_dphi(phi_vldr3_DM.data(), dphi_x.data(), dphi_y.data(), dphi_z.data(), fvl_thread); + } + if(isstress_) + { + phi_op.phi_dot_dphi_r(phi_vldr3_DM.data(), dphi_x.data(), dphi_y.data(), dphi_z.data(), 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; + } + } + } +} + + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.h new file mode 100644 index 0000000000..aec01a835b --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl.h @@ -0,0 +1,52 @@ +#pragma once + +#include +#include +#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 : public Gint +{ + public: + Gint_fvl( + const int nspin, + const std::vector& vr_eff, + const std::vector*>& DMR_vec, + const bool isforce, + const bool isstress, + ModuleBase::matrix* fvl, + ModuleBase::matrix* svl) + : nspin_(nspin), vr_eff_(vr_eff), 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 vr_eff_; + std::vector*> DMR_vec_; + const bool isforce_; + const bool isstress_; + + // output + ModuleBase::matrix* fvl_; + ModuleBase::matrix* svl_; + + // intermediate variables + std::vector>> DMRGint_vec_; + + double dr3_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.cpp new file mode 100644 index 0000000000..852669b264 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.cpp @@ -0,0 +1,132 @@ +#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(); + } +} + +void Gint_fvl_meta::cal_fvl_svl_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + std::vector phi_vldr3_DM; + std::vector dphi_x; + std::vector dphi_y; + std::vector dphi_z; + std::vector dphi_x_vldr3; + std::vector dphi_y_vldr3; + std::vector dphi_z_vldr3; + std::vector dphi_x_vldr3_DM; + std::vector dphi_y_vldr3_DM; + std::vector dphi_z_vldr3_DM; + std::vector ddphi_xx; + std::vector ddphi_xy; + std::vector ddphi_xz; + std::vector ddphi_yy; + std::vector ddphi_yz; + std::vector ddphi_zz; + 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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + phi_vldr3_DM.resize(phi_len); + dphi_x.resize(phi_len); + dphi_y.resize(phi_len); + dphi_z.resize(phi_len); + dphi_x_vldr3.resize(phi_len); + dphi_y_vldr3.resize(phi_len); + dphi_z_vldr3.resize(phi_len); + dphi_x_vldr3_DM.resize(phi_len); + dphi_y_vldr3_DM.resize(phi_len); + dphi_z_vldr3_DM.resize(phi_len); + ddphi_xx.resize(phi_len); + ddphi_xy.resize(phi_len); + ddphi_xz.resize(phi_len); + ddphi_yy.resize(phi_len); + ddphi_yz.resize(phi_len); + ddphi_zz.resize(phi_len); + phi_op.set_phi_dphi(phi.data(), dphi_x.data(), dphi_y.data(), dphi_z.data()); + phi_op.set_ddphi(ddphi_xx.data(), ddphi_xy.data(), ddphi_xz.data(), + ddphi_yy.data(), ddphi_yz.data(), ddphi_zz.data()); + for (int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_vldr3(vr_eff_[is], dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_x.data(), dphi_x_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_y.data(), dphi_y_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_z.data(), dphi_z_vldr3.data()); + phi_op.phi_mul_dm(phi_vldr3.data(), *DMRGint_vec_[is], false, phi_vldr3_DM.data()); + phi_op.phi_mul_dm(dphi_x_vldr3.data(), *DMRGint_vec_[is], false, dphi_x_vldr3_DM.data()); + phi_op.phi_mul_dm(dphi_y_vldr3.data(), *DMRGint_vec_[is], false, dphi_y_vldr3_DM.data()); + phi_op.phi_mul_dm(dphi_z_vldr3.data(), *DMRGint_vec_[is], false, dphi_z_vldr3_DM.data()); + if(isforce_) + { + phi_op.phi_dot_dphi(phi_vldr3_DM.data(), dphi_x.data(), dphi_y.data(), dphi_z.data(), fvl_thread); + phi_op.phi_dot_dphi(dphi_x_vldr3_DM.data(), ddphi_xx.data(), ddphi_xy.data(), ddphi_xz.data(), fvl_thread); + phi_op.phi_dot_dphi(dphi_y_vldr3_DM.data(), ddphi_xy.data(), ddphi_yy.data(), ddphi_yz.data(), fvl_thread); + phi_op.phi_dot_dphi(dphi_z_vldr3_DM.data(), ddphi_xz.data(), ddphi_yz.data(), ddphi_zz.data(), fvl_thread); + } + if(isstress_) + { + phi_op.phi_dot_dphi_r(phi_vldr3_DM.data(), dphi_x.data(), dphi_y.data(), dphi_z.data(), svl_thread); + phi_op.phi_dot_dphi_r(dphi_x_vldr3_DM.data(), ddphi_xx.data(), ddphi_xy.data(), ddphi_xz.data(), svl_thread); + phi_op.phi_dot_dphi_r(dphi_y_vldr3_DM.data(), ddphi_xy.data(), ddphi_yy.data(), ddphi_yz.data(), svl_thread); + phi_op.phi_dot_dphi_r(dphi_z_vldr3_DM.data(), ddphi_xz.data(), ddphi_yz.data(), ddphi_zz.data(), 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 \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.h new file mode 100644 index 0000000000..063e7971d1 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_fvl_meta.h @@ -0,0 +1,53 @@ +#pragma once + +#include +#include +#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& vr_eff, + const std::vector& vofk, + const std::vector*>& 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 vr_eff_; + std::vector vofk_; + std::vector*> DMR_vec_; + const bool isforce_; + const bool isstress_; + + // output + ModuleBase::matrix* fvl_; + ModuleBase::matrix* svl_; + + // intermediate variables + std::vector>> DMRGint_vec_; + + double dr3_; +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_helper.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_helper.h new file mode 100644 index 0000000000..b1f72c5a11 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_helper.h @@ -0,0 +1,63 @@ +#pragma once + +#include +#include +#include "gint_type.h" +#include "module_base/timer.h" + +template +std::shared_ptr toConstSharedPtr(std::shared_ptr ptr) { + return std::static_pointer_cast(ptr); +} + + +inline int index3Dto1D(const int id_x, const int id_y, const int id_z, + const int dim_x, const int dim_y, const int dim_z) +{ + return id_z + id_y * dim_z + id_x * dim_y * dim_z; +}; + +inline Vec3i index1Dto3D(const int index_1d, + const int dim_x, const int dim_y, const int dim_z) +{ + int id_x = index_1d / (dim_y * dim_z); + int id_y = (index_1d - id_x * dim_y * dim_z) / dim_z; + int id_z = index_1d % dim_z; + return Vec3i(id_x, id_y, id_z); +}; + +// if exponent is an integer between 0 and 5 (the most common cases in gint) and +// and exp is a variable that cannot be determined at compile time (which means the compiler cannot optimize the code), +// pow_int is much faster than std::pow +inline double pow_int(const double base, const int exp) +{ + switch (exp) + { + case 0: + return 1.0; + case 1: + return base; + case 2: + return base * base; + case 3: + return base * base * base; + case 4: + return base * base * base * base; + case 5: + return base * base * base * base * base; + default: + double result = std::pow(base, exp); + return result; + } +}; + +inline int floor_div(const int a, const int b) +{ + // a ^ b < 0 means a and b have different signs + return a / b - (a % b != 0 && (a ^ b) < 0); +}; + +inline int ceil_div(const int a, const int b) +{ + return a / b + (a % b != 0 && (a ^ b) > 0); +}; \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.cpp new file mode 100644 index 0000000000..1443ad05e9 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.cpp @@ -0,0 +1,212 @@ +#include +#include +#include "module_parameter/parameter.h" +#include "module_base/timer.h" +#include "gint_info.h" +#include "gint_type.h" + +namespace ModuleGint +{ + +GintInfo::GintInfo( + int nbx, int nby, int nbz, + int nmx, int nmy, int nmz, + int startidx_bx, int startidx_by, int startidx_bz, + int nbx_local, int nby_local, int nbz_local, + const Numerical_Orbital* Phi, + const UnitCell& ucell, Grid_Driver& gd) + : ucell_(&ucell) +{ + // initialize the unitcell information + unitcell_info_ = std::make_shared(ucell_->a1 * ucell_->lat0, ucell_->a2 * ucell_->lat0, ucell_->a3 * ucell_->lat0, + nbx, nby, nbz, nmx, nmy, nmz); + + biggrid_info_ = unitcell_info_->get_bgrid_info(); + meshgrid_info_ = biggrid_info_->get_mgrid_info(); + + // initialize the divide information + divide_info_ = std::make_shared(startidx_bx, startidx_by, startidx_bz, + nbx_local, nby_local, nbz_local, unitcell_info_, false); + + // initialize the localcell information + localcell_info_ = divide_info_->get_localcell_info(); + + // initialize the biggrids + BigGrid::init_localcell_info(localcell_info_); + BigGrid::init_unitcell_info(unitcell_info_); + BigGrid::init_bgrid_info(biggrid_info_); + + for (int i = 0; i < localcell_info_->get_bgrids_num(); i++) + { + biggrids_.push_back(std::make_shared(i)); + } + + // initialize the atoms + init_atoms_(ucell_->ntype, ucell_->atoms, Phi); + + // initialize the ijr_info + // this step needs to be done after init_atoms_, because it requires the information of is_atom_on_bgrid + init_ijr_info_(ucell, gd); +} + +template +std::shared_ptr> GintInfo::get_hr(int npol) const +{ + auto p_hr = std::make_shared>(ucell_->nat); + if(PARAM.inp.gamma_only) + { + p_hr->fix_gamma(); + } + p_hr->insert_ijrs(&ijr_info_, *ucell_, npol); + p_hr->allocate(nullptr, true); + return p_hr; +} + +void GintInfo::init_atoms_(int ntype, const Atom* atoms, const Numerical_Orbital* Phi) +{ + ModuleBase::timer::tick("GintInfo", "init_atoms"); + int iat = 0; + is_atom_in_proc_.resize(ucell_->nat, false); + atoms_.resize(ucell_->nat); + +// TODO: USE OPENMP TO PARALLELIZE THIS LOOP + for(int i = 0; i < ntype; i++) + { + const auto& atom = atoms[i]; + const auto *orb = &Phi[i]; + + // rcut extends to the maximum big grids in x, y, z directions + Vec3i ext_bgrid = biggrid_info_->max_ext_bgrid_num(atom.Rcut); + + for(int j = 0; j < atom.na; j++) + { + Vec3d fraction; + fraction.x = atom.taud[j].x * unitcell_info_->get_nbx(); + fraction.y = atom.taud[j].y * unitcell_info_->get_nby(); + fraction.z = atom.taud[j].z * unitcell_info_->get_nbz(); + const Vec3i atom_bgrid_idx(static_cast(fraction.x), + static_cast(fraction.y), + static_cast(fraction.z)); + const Vec3d delta(fraction.x - atom_bgrid_idx.x, + fraction.y - atom_bgrid_idx.y, + fraction.z - atom_bgrid_idx.z); + const Vec3d tau_in_biggrid = biggrid_info_->get_cartesian_coord(delta); + + const Vec3i ucell_idx_atom = unitcell_info_->get_unitcell_idx(atom_bgrid_idx); + auto& r_to_atom = atoms_[iat]; + + for(int bgrid_x = atom_bgrid_idx.x - ext_bgrid.x; bgrid_x <= atom_bgrid_idx.x + ext_bgrid.x; bgrid_x++) + { + for(int bgrid_y = atom_bgrid_idx.y - ext_bgrid.y; bgrid_y <= atom_bgrid_idx.y + ext_bgrid.y; bgrid_y++) + { + for(int bgrid_z = atom_bgrid_idx.z - ext_bgrid.z; bgrid_z <= atom_bgrid_idx.z + ext_bgrid.z; bgrid_z++) + { + // get the extended biggrid idx of the affected biggrid + const Vec3i ext_bgrid_idx(bgrid_x, bgrid_y, bgrid_z); + const Vec3i normal_bgrid_idx = unitcell_info_->map_ext_idx_to_ucell(ext_bgrid_idx); + if(localcell_info_->is_bgrid_in_lcell(normal_bgrid_idx) == false) + { + continue; + } + const int bgrid_local_idx = localcell_info_->get_bgrid_local_idx_1D(normal_bgrid_idx); + // get the unitcell idx of the big grid + const Vec3i ucell_idx_bgrid = unitcell_info_->get_unitcell_idx(ext_bgrid_idx); + + // The index of the unitcell containing the biggrid relative to the unitcell containing the atom. + const Vec3i ucell_idx_relative = ucell_idx_bgrid - ucell_idx_atom; + auto it = r_to_atom.find(ucell_idx_relative); + // if the gint_atom is not in the map, + // it means this is the first time we find this atom may affect some biggrids, + // add it to the r_to_atom map + if(it == r_to_atom.end()) + { + Vec3i ext_atom_bgrid_idx(atom_bgrid_idx.x - ucell_idx_bgrid.x * unitcell_info_->get_nbx(), + atom_bgrid_idx.y - ucell_idx_bgrid.y * unitcell_info_->get_nby(), + atom_bgrid_idx.z - ucell_idx_bgrid.z * unitcell_info_->get_nbz()); + r_to_atom.insert(std::make_pair(ucell_idx_relative, + GintAtom(&atom, j, iat, ext_atom_bgrid_idx, ucell_idx_relative, tau_in_biggrid, orb))); + } + if(biggrids_[bgrid_local_idx]->is_atom_on_bgrid(&r_to_atom.at(ucell_idx_relative))) + { + biggrids_[bgrid_local_idx]->add_atom(&r_to_atom.at(ucell_idx_relative)); + is_atom_in_proc_[iat] = true; + } + } + } + } + iat++; + } + } + ModuleBase::timer::tick("GintInfo", "init_atoms"); +} + +void GintInfo::init_ijr_info_(const UnitCell& ucell, Grid_Driver& gd) +{ + HContainer hRGint_local(ucell.nat); + // prepare the row_index and col_index for construct AtomPairs, they are + // same, name as orb_index + std::vector orb_index(ucell.nat + 1); + orb_index[0] = 0; + for (int i = 1; i < orb_index.size(); i++) { + int type = ucell.iat2it[i - 1]; + orb_index[i] = orb_index[i - 1] + ucell.atoms[type].nw; + } + + for (int T1 = 0; T1 < ucell.ntype; ++T1) { + const Atom* atom1 = &(ucell.atoms[T1]); + for (int I1 = 0; I1 < atom1->na; ++I1) { + auto& tau1 = atom1->tau[I1]; + const int iat1 = ucell.itia2iat(T1, I1); + // whether this atom is in this processor. + if (this->is_atom_in_proc_[iat1]) { + gd.Find_atom(ucell, tau1, T1, I1); + for (int ad = 0; ad < gd.getAdjacentNum() + 1; ++ad) { + const int T2 = gd.getType(ad); + const int I2 = gd.getNatom(ad); + const int iat2 = ucell.itia2iat(T2, I2); + const Atom* atom2 = &(ucell.atoms[T2]); + + // NOTE: hRGint wil save total number of atom pairs, + // if only upper triangle is saved, the lower triangle will + // be lost in 2D-block parallelization. if the adjacent atom + // is in this processor. + if (this->is_atom_in_proc_[iat2]) { + Vec3d dtau = gd.getAdjacentTau(ad) - tau1; + double distance = dtau.norm() * ucell.lat0; + double rcut = atom1->Rcut + atom2->Rcut; + + // if(distance < rcut) + // mohan reset this 2013-07-02 in Princeton + // we should make absolutely sure that the distance is + // smaller than rcuts[it] this should be consistant + // with LCAO_nnr::cal_nnrg function typical example : 7 + // Bohr cutoff Si orbital in 14 Bohr length of cell. + // distance = 7.0000000000000000 + // rcuts[it] = 7.0000000000000008 + if (distance < rcut - 1.0e-15) { + // calculate R index + auto& R_index = gd.getBox(ad); + // insert this atom-pair into this->hRGint + hamilt::AtomPair tmp_atom_pair( + iat1, + iat2, + R_index.x, + R_index.y, + R_index.z, + orb_index.data(), + orb_index.data(), + ucell.nat); + hRGint_local.insert_pair(tmp_atom_pair); + } + } + } + } + } + } + this->ijr_info_ = hRGint_local.get_ijr_info(); + return; +} + +template std::shared_ptr> GintInfo::get_hr(int npol) const; +template std::shared_ptr>> GintInfo::get_hr>(int npol) const; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.h new file mode 100644 index 0000000000..326f35a21b --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_info.h @@ -0,0 +1,82 @@ +#pragma once + +#include +#include +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" +#include "module_cell/atom_spec.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint_type.h" +#include "big_grid.h" +#include "gint_atom.h" +#include "unitcell_info.h" +#include "localcell_info.h" +#include "divide_info.h" + +namespace ModuleGint +{ + +class GintInfo +{ + public: + // constructor + GintInfo( + int nbx, int nby, int nbz, + int nmx, int nmy, int nmz, + int startidx_bx, int startidx_by, int startidx_bz, + int nbx_local, int nby_local, int nbz_local, + const Numerical_Orbital* Phi, + const UnitCell& ucell, Grid_Driver& gd); + + // getter functions + std::vector> get_biggrids() const { return biggrids_; }; + double get_local_mgrid_num() const { return localcell_info_->get_mgrids_num(); }; + double get_mgrid_volume() const { return meshgrid_info_->get_volume(); }; + + //========================================= + // functions about hcontainer + //========================================= + template + std::shared_ptr> get_hr(int npol = 1) const; + + private: + // initialize the atoms + void init_atoms_(int ntype, const Atom* atoms, const Numerical_Orbital* Phi); + + // initialize the ijr_info + void init_ijr_info_(const UnitCell& ucell, Grid_Driver& gd); + + const UnitCell* ucell_; + + // the unitcell information + std::shared_ptr unitcell_info_; + + // the biggrid information + std::shared_ptr biggrid_info_; + + // the meshgrid information + std::shared_ptr meshgrid_info_; + + // the divide information + std::shared_ptr divide_info_; + + // the localcell information + std::shared_ptr localcell_info_; + + // the big grids on this processor + std::vector> biggrids_; + + // the total atoms in the unitcell(include extended unitcell) on this processor + // atoms[iat][Vec3i] is the atom with index iat in the unitcell with index Vec3i + // Note: Since GintAtom does not implement a default constructor, + // the map should not be accessed using [], but rather using the at function + std::vector> atoms_; + + // if the iat-th(global index) atom is in this processor, return true + std::vector is_atom_in_proc_; + + // format for storing atomic pair information in hcontainer, used for initializing hcontainer + std::vector ijr_info_; +}; + +} // namespace ModuleGint diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.cpp new file mode 100644 index 0000000000..61f99f9118 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.cpp @@ -0,0 +1,110 @@ +#include "gint_interface.h" +#include "module_base/timer.h" +#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" + +namespace ModuleGint +{ + +void cal_gint_vl( + const double* vr_eff, + HContainer* hR) +{ + ModuleBase::timer::tick("Gint", "cal_gint_vl"); + Gint_vl gint_vl(vr_eff, hR); + gint_vl.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_vl"); +} + +void cal_gint_vl( + std::vector vr_eff, + HContainer>* hR) +{ + ModuleBase::timer::tick("Gint", "cal_gint_vl"); + Gint_vl_nspin4 gint_vl_nspin4(vr_eff, hR); + gint_vl_nspin4.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_vl"); +} + +void cal_gint_vl_metagga( + const double* vr_eff, + const double* vfork, + HContainer* hR) +{ + ModuleBase::timer::tick("Gint", "cal_gint_vl_metagga"); + Gint_vl_metagga gint_vl_metagga(vr_eff, vfork, hR); + gint_vl_metagga.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_vl_metagga"); +} + +void cal_gint_vl_metagga( + std::vector vr_eff, + std::vector vofk, + HContainer>* hR) +{ + ModuleBase::timer::tick("Gint", "cal_gint_vl_metagga"); + Gint_vl_metagga_nspin4 gint_vl_metagga_nspin4(vr_eff, vofk, hR); + gint_vl_metagga_nspin4.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_vl_metagga"); +} + +void cal_gint_rho( + const std::vector*>& DMR_vec, + const int nspin, + double **rho) +{ + ModuleBase::timer::tick("Gint", "cal_gint_rho"); + Gint_rho gint_rho(DMR_vec, nspin, rho); + gint_rho.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_rho"); +} + +void cal_gint_tau( + const std::vector*>& DMR_vec, + const int nspin, + double** tau) +{ + ModuleBase::timer::tick("Gint", "cal_gint_tau"); + Gint_tau gint_tau(DMR_vec, nspin, tau); + gint_tau.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_tau"); +} + +void cal_gint_fvl( + const int nspin, + const std::vector& vr_eff, + const std::vector*>& DMR_vec, + const bool isforce, + const bool isstress, + ModuleBase::matrix* fvl, + ModuleBase::matrix* svl) +{ + ModuleBase::timer::tick("Gint", "cal_gint_fvl"); + Gint_fvl gint_fvl(nspin, vr_eff, DMR_vec, isforce, isstress, fvl, svl); + gint_fvl.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_fvl"); +} + +void cal_gint_fvl_meta( + const int nspin, + const std::vector& vr_eff, + const std::vector& vofk, + const std::vector*>& DMR_vec, + const bool isforce, + const bool isstress, + ModuleBase::matrix* fvl, + ModuleBase::matrix* svl) +{ + ModuleBase::timer::tick("Gint", "cal_gint_fvl_meta"); + Gint_fvl_meta gint_fvl_meta(nspin, vr_eff, vofk, DMR_vec, isforce, isstress, fvl, svl); + gint_fvl_meta.cal_gint(); + ModuleBase::timer::tick("Gint", "cal_gint_fvl_meta"); +} + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.h new file mode 100644 index 0000000000..a4c33bb14a --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_interface.h @@ -0,0 +1,59 @@ +#pragma once +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint_type.h" + + +namespace ModuleGint +{ + +void cal_gint_vl( + const double* vr_eff, + HContainer* hR); + +void cal_gint_vl( + std::vector vr_eff, + HContainer>* hR); + +void cal_gint_vl_metagga( + const double* vr_eff, + const double* vfork, + HContainer* hR); + +void cal_gint_vl_metagga( + std::vector vr_eff, + std::vector vofk, + HContainer>* hR); + +void cal_gint_rho( + const std::vector*>& DMR_vec, + const int nspin, + double **rho); + +void cal_gint_tau( + const std::vector*>& DMR_vec, + const int nspin, + double**tau); + +void cal_gint_fvl( + const int nspin, + const std::vector& vr_eff, + const std::vector*>& DMR_vec, + const bool isforce, + const bool isstress, + ModuleBase::matrix* fvl, + ModuleBase::matrix* svl); + +void cal_gint_fvl_meta( + const int nspin, + const std::vector& vr_eff, + const std::vector& vofk, + const std::vector*>& DMR_vec, + const bool isforce, + const bool isstress, + ModuleBase::matrix* fvl, + ModuleBase::matrix* svl); + + + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.cpp new file mode 100644 index 0000000000..00655353bc --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.cpp @@ -0,0 +1,54 @@ +#include "module_base/global_function.h" +#include "gint_rho.h" +#include "gint_common.h" +#include "phi_operator.h" + +namespace ModuleGint +{ + +void Gint_rho::cal_gint() +{ + init_DMRGint_(); + transfer_DM_to_DMGint(gint_info_, DMR_vec_, DMRGint_vec_); + cal_rho_(); +} + +void Gint_rho::init_DMRGint_() +{ + DMRGint_vec_.resize(nspin_); + for (int is = 0; is < nspin_; is++) + { + DMRGint_vec_[is] = gint_info_->get_hr(); + } +} + +void Gint_rho::cal_rho_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_DMR; +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_DMR.resize(phi_len); + phi_op.set_phi(phi.data()); + for (int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_dm(phi.data(), *DMRGint_vec_[is], true, phi_DMR.data()); + phi_op.phi_dot_phi_dm(phi.data(), phi_DMR.data(), rho_[is]); + } + } + } +} + + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.h new file mode 100644 index 0000000000..227f604f84 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_rho.h @@ -0,0 +1,41 @@ +#pragma once + +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_rho : public Gint +{ + public: + Gint_rho( + const std::vector*>& DMR_vec, + const int nspin, + double **rho) + : DMR_vec_(DMR_vec), nspin_(nspin), rho_(rho) {}; + + void cal_gint() override; + + private: + void init_DMRGint_(); + + void cal_rho_(); + + // input + const std::vector*> DMR_vec_; + const int nspin_; + + // output + double **rho_; + + //======================== + // Intermediate variables + //======================== + std::vector>> DMRGint_vec_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.cpp new file mode 100644 index 0000000000..89374fbf00 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.cpp @@ -0,0 +1,65 @@ +#include "module_base/global_function.h" +#include "gint_tau.h" +#include "gint_common.h" +#include "phi_operator.h" + +namespace ModuleGint +{ + +void Gint_tau::cal_gint() +{ + init_DMRGint_(); + transfer_DM_to_DMGint(gint_info_, DMR_vec_, DMRGint_vec_); + cal_tau_(); +} + +void Gint_tau::init_DMRGint_() +{ + DMRGint_vec_.resize(nspin_); + for (int is = 0; is < nspin_; is++) + { + DMRGint_vec_[is] = gint_info_->get_hr(); + } +} + +void Gint_tau::cal_tau_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector dphi_x; + std::vector dphi_y; + std::vector dphi_z; + std::vector dphi_x_DM; + std::vector dphi_y_DM; + std::vector dphi_z_DM; +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + dphi_x.resize(phi_len); + dphi_y.resize(phi_len); + dphi_z.resize(phi_len); + dphi_x_DM.resize(phi_len); + dphi_y_DM.resize(phi_len); + dphi_z_DM.resize(phi_len); + phi_op.set_phi_dphi(nullptr, dphi_x.data(), dphi_y.data(), dphi_z.data()); + for (int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_dm(dphi_x.data(), *DMRGint_vec_[is], true, dphi_x_DM.data()); + phi_op.phi_mul_dm(dphi_y.data(), *DMRGint_vec_[is], true, dphi_y_DM.data()); + phi_op.phi_mul_dm(dphi_z.data(), *DMRGint_vec_[is], true, dphi_z_DM.data()); + phi_op.phi_dot_phi_dm(dphi_x.data(), dphi_x_DM.data(), kin_[is]); + phi_op.phi_dot_phi_dm(dphi_y.data(), dphi_y_DM.data(), kin_[is]); + phi_op.phi_dot_phi_dm(dphi_z.data(), dphi_z_DM.data(), kin_[is]); + } + } + } +} + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.h new file mode 100644 index 0000000000..65bb4620c7 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_tau.h @@ -0,0 +1,40 @@ +#pragma once +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_tau : public Gint +{ + public: + Gint_tau( + const std::vector*>& DMR_vec, + const int nspin, + double** tau) + : DMR_vec_(DMR_vec), nspin_(nspin), kin_(tau) {}; + + void cal_gint() override; + + private: + void init_DMRGint_(); + + void cal_tau_(); + + // input + const std::vector*> DMR_vec_; + const int nspin_; + + // output + double **kin_; + + //======================== + // Intermediate variables + //======================== + std::vector>> DMRGint_vec_; +}; + +} // namespace ModuleGint diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_type.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_type.h new file mode 100644 index 0000000000..c6e9369ad2 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_type.h @@ -0,0 +1,12 @@ +#pragma once + +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_base/vector3.h" +#include "module_base/matrix3.h" + +using Matrix3 = ModuleBase::Matrix3; +using Vec3d = ModuleBase::Vector3; +using Vec3i = ModuleBase::Vector3; + +template +using HContainer = hamilt::HContainer; \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.cpp new file mode 100644 index 0000000000..4553f73e76 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.cpp @@ -0,0 +1,60 @@ +#include "module_base/blas_connector.h" +#include "gint_common.h" +#include "gint_vl.h" +#include "phi_operator.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +void Gint_vl::cal_gint() +{ + init_hRGint_(); + cal_hRGint_(); + compose_hRGint(hRGint_); + transfer_hRGint_to_hR(toConstSharedPtr(hRGint_), hR_); +} + +//======================== +// Private functions +//======================== + +void Gint_vl::init_hRGint_() +{ + hRGint_ = gint_info_->get_hr(); +} + +void Gint_vl::cal_hRGint_() +{ +// be careful!! +// each thread will have a copy of hRGint_, this may cause a lot of memory usage +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + HContainer hRGint_local(*hRGint_); +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + phi_op.set_phi(phi.data()); + phi_op.phi_mul_vldr3(vr_eff_, dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_phi_vldr3(phi.data(), phi_vldr3.data(), &hRGint_local); + } +#pragma omp critical + { + BlasConnector::axpy(hRGint_local.get_nnr(), 1.0, hRGint_local.get_wrapper(), + 1, hRGint_->get_wrapper(), 1); + } + } +} + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.h new file mode 100644 index 0000000000..b68cb1a4a7 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl.h @@ -0,0 +1,44 @@ +#pragma once + +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_vl : public Gint +{ + public: + Gint_vl( + const double* vr_eff, + HContainer* hR) + : vr_eff_(vr_eff), hR_(hR), dr3_(gint_info_->get_mgrid_volume()){}; + + void cal_gint() override; + + private: + + void init_hRGint_(); + + // note that only the upper triangle matrix of hR is calculated + // that's why we need compose_hRGint() to fill the lower triangle matrix. + void cal_hRGint_(); + + // input + const double* vr_eff_; + + // output + HContainer* hR_; + + //======================== + // Intermediate variables + //======================== + double dr3_; + + std::shared_ptr> hRGint_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.cpp new file mode 100644 index 0000000000..2c5351c338 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.cpp @@ -0,0 +1,78 @@ +#include "module_base/blas_connector.h" +#include "gint_common.h" +#include "gint_vl_metagga.h" +#include "phi_operator.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +void Gint_vl_metagga::cal_gint() +{ + init_hRGint_(); + cal_hRGint_(); + compose_hRGint(hRGint_); + transfer_hRGint_to_hR(toConstSharedPtr(hRGint_), hR_); +} + +//======================== +// Private functions +//======================== + +void Gint_vl_metagga::init_hRGint_() +{ + hRGint_ = gint_info_->get_hr(); +} + +void Gint_vl_metagga::cal_hRGint_() +{ +// be careful!! +// each thread will have a copy of hRGint_, this may cause a lot of memory usage +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + std::vector dphi_x; + std::vector dphi_y; + std::vector dphi_z; + std::vector dphi_x_vldr3; + std::vector dphi_y_vldr3; + std::vector dphi_z_vldr3; + HContainer hRGint_local(*hRGint_); +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + dphi_x.resize(phi_len); + dphi_y.resize(phi_len); + dphi_z.resize(phi_len); + dphi_x_vldr3.resize(phi_len); + dphi_y_vldr3.resize(phi_len); + dphi_z_vldr3.resize(phi_len); + phi_op.set_phi_dphi(phi.data(), dphi_x.data(), dphi_y.data(), dphi_z.data()); + phi_op.phi_mul_vldr3(vr_eff_, dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_, dr3_, dphi_x.data(), dphi_x_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_, dr3_, dphi_y.data(), dphi_y_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_, dr3_, dphi_z.data(), dphi_z_vldr3.data()); + phi_op.phi_mul_phi_vldr3(phi.data(), phi_vldr3.data(), &hRGint_local); + phi_op.phi_mul_phi_vldr3(dphi_x.data(), dphi_x_vldr3.data(), &hRGint_local); + phi_op.phi_mul_phi_vldr3(dphi_y.data(), dphi_y_vldr3.data(), &hRGint_local); + phi_op.phi_mul_phi_vldr3(dphi_z.data(), dphi_z_vldr3.data(), &hRGint_local); + } +#pragma omp critical + { + BlasConnector::axpy(hRGint_local.get_nnr(), 1.0, hRGint_local.get_wrapper(), + 1, hRGint_->get_wrapper(), 1); + } + } +} + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.h new file mode 100644 index 0000000000..4fef9c85af --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga.h @@ -0,0 +1,47 @@ +#pragma once + +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_vl_metagga : public Gint +{ + public: + Gint_vl_metagga( + const double* vr_eff, + const double* vofk, + HContainer* hR) + : vr_eff_(vr_eff), vofk_(vofk), hR_(hR), dr3_(gint_info_->get_mgrid_volume()){}; + + void cal_gint() override; + + private: + + void init_hRGint_(); + + // note that only the upper triangle matrix of hR is calculated + // that's why we need compose_hRGint() to fill the lower triangle matrix. + void cal_hRGint_(); + + // input + const double* vr_eff_; + const double* vofk_; + + // output + HContainer* hR_; + + //======================== + // Intermediate variables + //======================== + double dr3_; + + std::shared_ptr> hRGint_; + +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.cpp new file mode 100644 index 0000000000..9ee94b466c --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.cpp @@ -0,0 +1,88 @@ +#include "module_base/global_function.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_base/blas_connector.h" +#include "gint_common.h" +#include "gint_vl_metagga_nspin4.h" +#include "phi_operator.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +void Gint_vl_metagga_nspin4::cal_gint() +{ + init_hRGint_(); + cal_hRGint_(); + compose_hRGint(hRGint_part_, hRGint_full_); + transfer_hRGint_to_hR(toConstSharedPtr(hRGint_full_), hR_); +} + +void Gint_vl_metagga_nspin4::init_hRGint_() +{ + hRGint_part_.resize(nspin_); + for(int i = 0; i < nspin_; i++) + { + hRGint_part_[i] = gint_info_->get_hr(); + } + const int npol = 2; + hRGint_full_ = gint_info_->get_hr>(npol); +} + +void Gint_vl_metagga_nspin4::cal_hRGint_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + std::vector dphi_x; + std::vector dphi_y; + std::vector dphi_z; + std::vector dphi_x_vldr3; + std::vector dphi_y_vldr3; + std::vector dphi_z_vldr3; + std::vector> hRGint_part_thread(nspin_, *hRGint_part_[0]); +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + dphi_x.resize(phi_len); + dphi_y.resize(phi_len); + dphi_z.resize(phi_len); + dphi_x_vldr3.resize(phi_len); + dphi_y_vldr3.resize(phi_len); + dphi_z_vldr3.resize(phi_len); + phi_op.set_phi_dphi(phi.data(), dphi_x.data(), dphi_y.data(), dphi_z.data()); + for(int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_vldr3(vr_eff_[is], dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_x.data(), dphi_x_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_y.data(), dphi_y_vldr3.data()); + phi_op.phi_mul_vldr3(vofk_[is], dr3_, dphi_z.data(), dphi_z_vldr3.data()); + phi_op.phi_mul_phi_vldr3(phi.data(), phi_vldr3.data(), &hRGint_part_thread[is]); + phi_op.phi_mul_phi_vldr3(dphi_x.data(), dphi_x_vldr3.data(), &hRGint_part_thread[is]); + phi_op.phi_mul_phi_vldr3(dphi_y.data(), dphi_y_vldr3.data(), &hRGint_part_thread[is]); + phi_op.phi_mul_phi_vldr3(dphi_z.data(), dphi_z_vldr3.data(), &hRGint_part_thread[is]); + } + } +#pragma omp critical + { + for(int is = 0; is < nspin_; is++) + { + { + BlasConnector::axpy(hRGint_part_thread[is].get_nnr(), 1.0, hRGint_part_thread[is].get_wrapper(), + 1, hRGint_part_[is]->get_wrapper(), 1); + } + } + } + } +} + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.h new file mode 100644 index 0000000000..d8033e4be6 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_metagga_nspin4.h @@ -0,0 +1,45 @@ +#pragma once + +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_vl_metagga_nspin4 : public Gint +{ + public: + Gint_vl_metagga_nspin4( + std::vector vr_eff, + std::vector vofk, + HContainer>* hR) + : vr_eff_(vr_eff), vofk_(vofk), hR_(hR), dr3_(gint_info_->get_mgrid_volume()){}; + + void cal_gint() override; + + private: + void init_hRGint_(); + + void cal_hRGint_(); + + // input + std::vector vr_eff_; + std::vector vofk_; + // output + HContainer>* hR_; + + //======================== + // Intermediate variables + //======================== + const double dr3_; + + const int nspin_ = 4; + + std::vector>> hRGint_part_; + std::shared_ptr>> hRGint_full_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.cpp new file mode 100644 index 0000000000..b9501940bf --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.cpp @@ -0,0 +1,69 @@ +#include "module_base/global_function.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_base/blas_connector.h" +#include "gint_common.h" +#include "gint_vl_nspin4.h" +#include "phi_operator.h" +#include "gint_helper.h" + +namespace ModuleGint +{ +void Gint_vl_nspin4::cal_gint() +{ + init_hRGint_(); + cal_hRGint_(); + compose_hRGint(hRGint_part_, hRGint_full_); + transfer_hRGint_to_hR(toConstSharedPtr(hRGint_full_), hR_); +} + +void Gint_vl_nspin4::init_hRGint_() +{ + hRGint_part_.resize(nspin_); + for(int i = 0; i < nspin_; i++) + { + hRGint_part_[i] = gint_info_->get_hr(); + } + const int npol = 2; + hRGint_full_ = gint_info_->get_hr>(npol); +} + +void Gint_vl_nspin4::cal_hRGint_() +{ +#pragma omp parallel + { + PhiOperator phi_op; + std::vector phi; + std::vector phi_vldr3; + std::vector> hRGint_part_thread(nspin_, *hRGint_part_[0]); +#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); + const int phi_len = phi_op.get_rows() * phi_op.get_cols(); + phi.resize(phi_len); + phi_vldr3.resize(phi_len); + phi_op.set_phi(phi.data()); + for(int is = 0; is < nspin_; is++) + { + phi_op.phi_mul_vldr3(vr_eff_[is], dr3_, phi.data(), phi_vldr3.data()); + phi_op.phi_mul_phi_vldr3(phi.data(), phi_vldr3.data(), &hRGint_part_thread[is]); + } + } +#pragma omp critical + { + for(int is = 0; is < nspin_; is++) + { + { + BlasConnector::axpy(hRGint_part_thread[is].get_nnr(), 1.0, hRGint_part_thread[is].get_wrapper(), + 1, hRGint_part_[is]->get_wrapper(), 1); + } + } + } + } +} + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.h b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.h new file mode 100644 index 0000000000..3f4a725459 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/gint_vl_nspin4.h @@ -0,0 +1,48 @@ +#pragma once + +#include +#include +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "gint.h" +#include "gint_info.h" + +namespace ModuleGint +{ + +class Gint_vl_nspin4 : public Gint +{ + public: + Gint_vl_nspin4( + std::vector vr_eff, + HContainer>* hR) + : vr_eff_(vr_eff), hR_(hR), dr3_(gint_info_->get_mgrid_volume()){}; + + void cal_gint() override; + + private: + + void init_hRGint_(); + + // note that only the upper triangle matrix of hR is calculated + // that's why we need compose_hRGint() to fill the lower triangle matrix. + void cal_hRGint_(); + + // input + std::vector vr_eff_; + + // output + HContainer>* hR_; + + //======================== + // Intermediate variables + //======================== + const double dr3_; + + const int nspin_ = 4; + + std::vector>> hRGint_part_; + std::shared_ptr>> hRGint_full_; + +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.cpp new file mode 100644 index 0000000000..da9b5ec3a3 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.cpp @@ -0,0 +1,119 @@ +#include "localcell_info.h" + +namespace ModuleGint +{ + LocalCellInfo::LocalCellInfo( + int startidx_bx, int startidx_by, int startidx_bz, + int nbx, int nby, int nbz, + std::shared_ptr unitcell_info) + : startidx_bx_(startidx_bx), startidx_by_(startidx_by), startidx_bz_(startidx_bz), + nbx_(nbx), nby_(nby), nbz_(nbz), nbxyz_(nbx*nby*nbz), + unitcell_info_(unitcell_info), biggrid_info_(unitcell_info->get_bgrid_info()) + { + startidx_mx_ = startidx_bx_ * biggrid_info_->get_nmx(); + startidx_my_ = startidx_by_ * biggrid_info_->get_nmy(); + startidx_mz_ = startidx_bz_ * biggrid_info_->get_nmz(); + nmx_ = nbx_ * biggrid_info_->get_nmx(); + nmy_ = nby_ * biggrid_info_->get_nmy(); + nmz_ = nbz_ * biggrid_info_->get_nmz(); + nmxyz_ = nmx_ * nmy_ * nmz_; + } + + //==================================================================== + // functions related to the big grid + //==================================================================== + + int LocalCellInfo::bgrid_idx_3Dto1D(const Vec3i index_3d) const + { + return index3Dto1D(index_3d.x, index_3d.y, index_3d.z, nbx_, nby_, nbz_); + } + + Vec3i LocalCellInfo::bgrid_idx_1Dto3D(const int index_1d) const + { + return index1Dto3D(index_1d, nbx_, nby_, nbz_); + } + + Vec3i LocalCellInfo::get_bgrid_global_idx_3D(const Vec3i index_3d) const + { + return Vec3i( + startidx_bx_ + index_3d.x, + startidx_by_ + index_3d.y, + startidx_bz_ + index_3d.z); + } + + Vec3i LocalCellInfo::get_bgrid_global_idx_3D(const int index_1d) const + { + return get_bgrid_global_idx_3D(bgrid_idx_1Dto3D(index_1d)); + } + + int LocalCellInfo::get_bgrid_global_idx_1D(const int index_1d) const + { + Vec3i ucell_idx_3d = get_bgrid_global_idx_3D(bgrid_idx_1Dto3D(index_1d)); + return unitcell_info_->bgrid_idx_3Dto1D(ucell_idx_3d); + } + + + Vec3i LocalCellInfo::get_bgrid_local_idx_3D(const Vec3i index_3d) const + { + int x = index_3d.x - startidx_bx_; + int y = index_3d.y - startidx_by_; + int z = index_3d.z - startidx_bz_; + assert(x >= 0 && x < nbx_); + assert(y >= 0 && y < nby_); + assert(z >= 0 && z < nbz_); + return Vec3i(x, y, z); + } + + int LocalCellInfo::get_bgrid_local_idx_1D(const Vec3i index_3d) const + { + return bgrid_idx_3Dto1D(get_bgrid_local_idx_3D(index_3d)); + } + + int LocalCellInfo::get_bgrid_local_idx_1D(const int index_1d) const + { + Vec3i idx_3d = unitcell_info_->bgrid_idx_1Dto3D(index_1d); + return bgrid_idx_3Dto1D(get_bgrid_local_idx_3D(idx_3d)); + } + + Vec3d LocalCellInfo::get_bgrid_global_coord_3D(const int index_1d) const + { + Vec3i ucell_idx_3d = get_bgrid_global_idx_3D(index_1d); + return unitcell_info_->get_bgrid_coord(ucell_idx_3d); + } + + bool LocalCellInfo::is_bgrid_in_lcell(const Vec3i index_3d) const + { + return (index_3d.x >= startidx_bx_ && index_3d.x < startidx_bx_ + nbx_ && + index_3d.y >= startidx_by_ && index_3d.y < startidx_by_ + nby_ && + index_3d.z >= startidx_bz_ && index_3d.z < startidx_bz_ + nbz_); + } + + //==================================================================== + // functions related to the meshgrid + //==================================================================== + + int LocalCellInfo::mgrid_idx_3Dto1D(const Vec3i index_3d) const + { + return index3Dto1D(index_3d.x, index_3d.y, index_3d.z, nmx_, nmy_, nmz_); + } + + Vec3i LocalCellInfo::mgrid_idx_1Dto3D(const int index_1d) const + { + return index1Dto3D(index_1d, nmx_, nmy_, nmz_); + } + + Vec3i LocalCellInfo::get_mgrid_global_idx_3D(const Vec3i index_3d) const + { + return Vec3i( + startidx_mx_ + index_3d.x, + startidx_my_ + index_3d.y, + startidx_mz_ + index_3d.z); + } + + int LocalCellInfo::get_mgrid_global_idx_1D(const int index_1d) const + { + Vec3i ucell_idx_3d = get_mgrid_global_idx_3D(mgrid_idx_1Dto3D(index_1d)); + return unitcell_info_->mgrid_idx_3Dto1D(ucell_idx_3d); + } + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h new file mode 100644 index 0000000000..a39042e558 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/localcell_info.h @@ -0,0 +1,126 @@ +#pragma once + +#include +#include "gint_type.h" +#include "unitcell_info.h" + +namespace ModuleGint +{ + +class LocalCellInfo +{ + public: + // constructor + LocalCellInfo( + int startidx_x, int startidx_y, int startidx_z, + int nbx, int nby, int nbz, + std::shared_ptr unitcell_info); + + // getter functions + const int get_startidx_bx() const { return startidx_bx_; }; + const int get_startidx_by() const { return startidx_by_; }; + const int get_startidx_bz() const { return startidx_bz_; }; + const int get_nbx() const { return nbx_; }; + const int get_nby() const { return nby_; }; + const int get_nbz() const { return nbz_; }; + const int get_bgrids_num() const { return nbxyz_; }; + const int get_mgrids_num() const { return nmxyz_; }; + std::shared_ptr get_unitcell_info() const { return unitcell_info_; }; + std::shared_ptr get_bgrid_info() const { return unitcell_info_->get_bgrid_info(); }; + + //==================================================================== + // functions related to the big grid + //==================================================================== + + // transform the 3D index of a big grid in the local cell to the 3D index in the local cell + int bgrid_idx_3Dto1D(const Vec3i index_3d) const; + + // transform the 1D index of a big grid in the local cell to the 1D index in the local cell + Vec3i bgrid_idx_1Dto3D(const int index_1d) const; + + // transform the 3D index of a big grid in the local cell to the 3D index in the unit cell + Vec3i get_bgrid_global_idx_3D(const Vec3i index_3d) const; + + // transform the 1D index of a big grid in the local cell to the 3D index in the unit cell + Vec3i get_bgrid_global_idx_3D(const int index_1d) const; + + // transform the 1D index of a big grid in the local cell to the 1D index in the unit cell + int get_bgrid_global_idx_1D(const int index_1d) const; + + // transform the 3D index of a big grid in the unit cell to the 3D index in the local cell + Vec3i get_bgrid_local_idx_3D(const Vec3i index_3d) const; + + // transform the 1D index of a big grid in the unit cell to the 1D index in the local cell + int get_bgrid_local_idx_1D(const int index_1d) const; + + // transform the 3D index of a big grid in the unit cell to the 1D index in the local cell + int get_bgrid_local_idx_1D(const Vec3i index_3d) const; + + // get the cartesian coordinate of a big grid in the unit cell from the 1D index + Vec3d get_bgrid_global_coord_3D(const int index_1d) const; + + // the input is the 3D index of a big grid in the unitcell + // return true if the big grid is in the local cell + bool is_bgrid_in_lcell(const Vec3i index_3d) const; + + + //==================================================================== + // functions related to the meshgrid + //==================================================================== + + // transform the 3D index of a meshgrid in the local cell to the 3D index in the local cell + int mgrid_idx_3Dto1D(const Vec3i index_3d) const; + + // transform the 1D index of a meshgrid in the local cell to the 1D index in the local cell + Vec3i mgrid_idx_1Dto3D(const int index_1d) const; + + // transform the 3D index of a meshgrid in the local cell to the 3D index in the unit cell + Vec3i get_mgrid_global_idx_3D(const Vec3i index_3d) const; + + // transform the 1D index of a meshgrid in the local cell to the 1D index in the unit cell + int get_mgrid_global_idx_1D(const int index_1d) const; + + private: + //==================================================================== + // information about the big grid + //==================================================================== + + // 3D index of the first big grid in the local cell within the unit cell + int startidx_bx_; + int startidx_by_; + int startidx_bz_; + + // Number of big grids in the local cell along the three basis vectors of the local cell + int nbx_; + int nby_; + int nbz_; + + // Total number of big grids in the local cell + int nbxyz_; + + //==================================================================== + // information about the meshgrid + //==================================================================== + + // 3D index of the first meshgrid in the local cell within the unit cell + int startidx_mx_; + int startidx_my_; + int startidx_mz_; + + // Number of meshgrids in the local cell along the three basis vectors of the local cell + int nmx_; + int nmy_; + int nmz_; + + // Total number of meshgrids in the local cell + int nmxyz_; + + // information about the Unitcell + std::shared_ptr unitcell_info_; + + // information about the big grid + std::shared_ptr biggrid_info_; + +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/meshgrid_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/meshgrid_info.h new file mode 100644 index 0000000000..581c9944f7 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/meshgrid_info.h @@ -0,0 +1,64 @@ +#pragma once + +#include "gint_type.h" +#include "module_cell/unitcell.h" + +namespace ModuleGint +{ + +class MeshGridInfo +{ + public: + // constructor + MeshGridInfo( + Vec3d meshgrid_vec1, + Vec3d meshgrid_vec2, + Vec3d meshgrid_vec3) + : meshgrid_vec1_(meshgrid_vec1), + meshgrid_vec2_(meshgrid_vec2), + meshgrid_vec3_(meshgrid_vec3) + { + // initialize the meshgrid_latvec0_ + meshgrid_latvec0_.e11 = meshgrid_vec1_.x; + meshgrid_latvec0_.e12 = meshgrid_vec1_.y; + meshgrid_latvec0_.e13 = meshgrid_vec1_.z; + + meshgrid_latvec0_.e21 = meshgrid_vec2_.x; + meshgrid_latvec0_.e22 = meshgrid_vec2_.y; + meshgrid_latvec0_.e23 = meshgrid_vec2_.z; + + meshgrid_latvec0_.e31 = meshgrid_vec3_.x; + meshgrid_latvec0_.e32 = meshgrid_vec3_.y; + meshgrid_latvec0_.e33 = meshgrid_vec3_.z; + + // initialize the GT matrix + meshgrid_GT_ = meshgrid_latvec0_.Inverse(); + + meshgrid_volume_ = std::abs(meshgrid_latvec0_.Det()); + }; + + double get_volume() const { return meshgrid_volume_; }; + Vec3d get_cartesian_coord(const Vec3i& index_3d) const { return index_3d * meshgrid_latvec0_; }; + Vec3d get_direct_coord(const Vec3d& cart_coord) const { return cart_coord * meshgrid_GT_; } + + private: + // basis vectors of meshgrid + Vec3d meshgrid_vec1_; + Vec3d meshgrid_vec2_; + Vec3d meshgrid_vec3_; + + // used to convert the (i, j, k) index of the meshgrid to the Cartesian coordinate + // if meshrid_vec1_ is row vector, + // then meshgrid_latvec0_ = [meshgrid_vec1_; meshgrid_vec2_; meshgrid_vec3_], + // (i, j, k) * meshgrid_latvec0_ = (x, y, z) + Matrix3 meshgrid_latvec0_; + + // used to convert the Cartesian coordinate to the (i, j, k) index of the mesh grid + // meshgrid_GT_ = meshgrid_latvec0_.Inverse() + // (x, y, z) * meshgrid_GT_ = (i, j, k) + Matrix3 meshgrid_GT_; + + double meshgrid_volume_; +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.cpp new file mode 100644 index 0000000000..6e33e3a903 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.cpp @@ -0,0 +1,317 @@ +#include "phi_operator.h" +#include "module_base/blas_connector.h" +#include "module_base/global_function.h" +#include "module_base/matrix.h" + +namespace ModuleGint +{ + +void PhiOperator::set_bgrid(std::shared_ptr biggrid) +{ + biggrid_ = biggrid; + rows_ = biggrid_->get_mgrids_num(); + cols_ = biggrid_->get_mgrid_phi_len(); + + biggrid_->set_atoms_startidx(atoms_startidx_); + biggrid_->set_atoms_phi_len(atoms_phi_len_); + biggrid_->set_mgrids_local_idx(meshgrids_local_idx_); + + // init is_atom_on_mgrid_ and atoms_relative_coords_ + int atoms_num = biggrid_->get_atoms_num(); + atoms_relative_coords_.resize(atoms_num); + is_atom_on_mgrid_.resize(atoms_num); + for(int i = 0; i < atoms_num; ++i) + { + biggrid_->set_atom_relative_coords(biggrid_->get_atom(i), atoms_relative_coords_[i]); + is_atom_on_mgrid_[i].resize(rows_); + for(int j = 0; j < rows_; ++j) + { + is_atom_on_mgrid_[i][j] = atoms_relative_coords_[i][j].norm() <= biggrid_->get_atom(i)->get_rcut(); + } + } + + // init atom_pair_start_end_idx_ + init_atom_pair_start_end_idx_(); +} + +void PhiOperator::set_phi(double* phi) const +{ + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const auto atom = biggrid_->get_atom(i); + atom->set_phi(atoms_relative_coords_[i], cols_, phi); + phi += atom->get_nw(); + } +} + +void PhiOperator::set_phi_dphi(double* phi, double* dphi_x, double* dphi_y, double* dphi_z) const +{ + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const auto atom = biggrid_->get_atom(i); + atom->set_phi_dphi(atoms_relative_coords_[i], cols_, phi, dphi_x, dphi_y, dphi_z); + if(phi != nullptr) + { + phi += atom->get_nw(); + } + dphi_x += atom->get_nw(); + dphi_y += atom->get_nw(); + dphi_z += atom->get_nw(); + } +} + +void PhiOperator::set_ddphi( + double* ddphi_xx, double* ddphi_xy, double* ddphi_xz, + double* ddphi_yy, double* ddphi_yz, double* ddphi_zz) const +{ + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const auto atom = biggrid_->get_atom(i); + atom->set_ddphi(atoms_relative_coords_[i], cols_, ddphi_xx, ddphi_xy, ddphi_xz, ddphi_yy, ddphi_yz, ddphi_zz); + ddphi_xx += atom->get_nw(); + ddphi_xy += atom->get_nw(); + ddphi_xz += atom->get_nw(); + ddphi_yy += atom->get_nw(); + ddphi_yz += atom->get_nw(); + ddphi_zz += atom->get_nw(); + } +} + +void PhiOperator::phi_mul_dm( + const double* phi, + const HContainer& DM, + const bool is_symm, double* phi_dm) const +{ + ModuleBase::GlobalFunc::ZEROS(phi_dm, rows_ * cols_); + // parameters for lapack subroutines + constexpr char side = 'L'; + constexpr char uplo = 'U'; + const char trans = 'N'; + const double alpha = 1.0; + const double beta = 1.0; + const double alpha1 = is_symm ? 2.0 : 1.0; + + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const auto atom_i = biggrid_->get_atom(i); + const auto r_i = atom_i->get_R(); + + if(is_symm) + { + const auto dm_mat = DM.find_matrix(atom_i->get_iat(), atom_i->get_iat(), 0, 0, 0); + dsymm_(&side, &uplo, &atoms_phi_len_[i], &rows_, &alpha, dm_mat->get_pointer(), &atoms_phi_len_[i], + &phi[0 * cols_ + atoms_startidx_[i]], &cols_, &beta, &phi_dm[0 * cols_ + atoms_startidx_[i]], &cols_); + } + + const int start = is_symm ? i + 1 : 0; + + for(int j = start; j < biggrid_->get_atoms_num(); ++j) + { + const auto atom_j = biggrid_->get_atom(j); + const auto r_j = atom_j->get_R(); + // FIXME may be r = r_j - r_i + const auto dm_mat = DM.find_matrix(atom_i->get_iat(), atom_j->get_iat(), r_i-r_j); + + // if dm_mat is nullptr, it means this atom pair does not affect any meshgrid in the unitcell + if(dm_mat == nullptr) + { + continue; + } + + int start_idx = get_atom_pair_start_end_idx_(i, j).first; + int end_idx = get_atom_pair_start_end_idx_(i, j).second; + const int len = end_idx - start_idx + 1; + + // if len<=0, it means this atom pair does not affect any meshgrid in this biggrid + if(len <= 0) + { + continue; + } + + dgemm_(&trans, &trans, &atoms_phi_len_[j], &len, &atoms_phi_len_[i], &alpha1, dm_mat->get_pointer(), &atoms_phi_len_[j], + &phi[start_idx * cols_ + atoms_startidx_[i]], &cols_, &beta, &phi_dm[start_idx * cols_ + atoms_startidx_[j]], &cols_); + } + } +} + +void PhiOperator::phi_mul_vldr3(const double* vl, const double dr3, const double* phi, double* result) const +{ + int idx = 0; + for(int i = 0; i < biggrid_->get_mgrids_num(); i++) + { + double vldr3_mgrid = vl[meshgrids_local_idx_[i]] * dr3; + for(int j = 0; j < cols_; j++) + { + result[idx] = phi[idx] * vldr3_mgrid; + idx++; + } + } +} + +void PhiOperator::phi_mul_phi_vldr3( + const double* phi, + const double* phi_vldr3, + HContainer* hr) const +{ + const char transa='N', transb='T'; + const double alpha=1, beta=1; + + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const auto atom_i = biggrid_->get_atom(i); + const auto& r_i = atom_i->get_R(); + const int iat_i = atom_i->get_iat(); + + for(int j = 0; j < biggrid_->get_atoms_num(); ++j) + { + const auto atom_j = biggrid_->get_atom(j); + const auto& r_j = atom_j->get_R(); + const int iat_j = atom_j->get_iat(); + + // only calculate the upper triangle matrix + if(iat_i > iat_j) + { + continue; + } + + // FIXME may be r = r_j - r_i + const auto result = hr->find_matrix(iat_i, iat_j, r_i-r_j); + + if(result == nullptr) + { + continue; + } + + int start_idx = get_atom_pair_start_end_idx_(i, j).first; + int end_idx = get_atom_pair_start_end_idx_(i, j).second; + const int len = end_idx - start_idx + 1; + + if(len <= 0) + { + continue; + } + + dgemm_(&transa, &transb, &atoms_phi_len_[j], &atoms_phi_len_[i], &len, &alpha, &phi_vldr3[start_idx * cols_ + atoms_startidx_[j]], + &cols_,&phi[start_idx * cols_ + atoms_startidx_[i]], &cols_, &beta, result->get_pointer(), &atoms_phi_len_[j]); + } + } +} + +void PhiOperator::phi_dot_phi_dm( + const double* phi, + const double* phi_dm, + double* rho) const +{ + const int inc = 1; + for(int i = 0; i < biggrid_->get_mgrids_num(); ++i) + { + rho[meshgrids_local_idx_[i]] += ddot_(&cols_, &phi[i * cols_], &inc, &phi_dm[i * cols_], &inc); + } +} + +void PhiOperator::phi_dot_dphi( + const double* phi, + const double* dphi_x, + const double* dphi_y, + const double* dphi_z, + ModuleBase::matrix *fvl) const +{ + for(int i = 0; i < biggrid_->get_atoms_num(); ++i) + { + const int iat = biggrid_->get_atom(i)->get_iat(); + const int start_idx = atoms_startidx_[i]; + const int phi_len = atoms_phi_len_[i]; + double rx = 0, ry = 0, rz = 0; + for(int j = 0; j < biggrid_->get_mgrids_num(); ++j) + { + for(int k = 0; k < phi_len; ++k) + { + int idx = j * cols_ + start_idx + k; + const double phi_val = phi[idx]; + rx += phi_val * dphi_x[idx]; + ry += phi_val * dphi_y[idx]; + rz += phi_val * dphi_z[idx]; + } + } + fvl[0](iat, 0) += rx * 2; + fvl[0](iat, 1) += ry * 2; + fvl[0](iat, 2) += rz * 2; + } +} + +void PhiOperator::phi_dot_dphi_r( + const double* phi, + const double* dphi_x, + const double* dphi_y, + const double* dphi_z, + ModuleBase::matrix *svl) const +{ + double sxx = 0, sxy = 0, sxz = 0, syy = 0, syz = 0, szz = 0; + for(int i = 0; i < biggrid_->get_mgrids_num(); ++i) + { + for(int j = 0; j < biggrid_->get_atoms_num(); ++j) + { + const int start_idx = atoms_startidx_[j]; + for(int k = 0; k < atoms_phi_len_[j]; ++k) + { + const int idx = i * cols_ + start_idx + k; + const Vec3d& r3 = atoms_relative_coords_[j][i]; + const double phi_val = phi[idx]; + sxx += phi_val * dphi_x[idx] * r3[0]; + sxy += phi_val * dphi_x[idx] * r3[1]; + sxz += phi_val * dphi_x[idx] * r3[2]; + syy += phi_val * dphi_y[idx] * r3[1]; + syz += phi_val * dphi_y[idx] * r3[2]; + szz += phi_val * dphi_z[idx] * r3[2]; + } + } + } + svl[0](0, 0) += sxx * 2; + svl[0](0, 1) += sxy * 2; + svl[0](0, 2) += sxz * 2; + svl[0](1, 1) += syy * 2; + svl[0](1, 2) += syz * 2; + svl[0](2, 2) += szz * 2; +} + + +//=============================== +// private methods +//=============================== +void PhiOperator::init_atom_pair_start_end_idx_() +{ + int atoms_num = biggrid_->get_atoms_num(); + atom_pair_start_end_idx_.resize(atoms_num * (atoms_num + 1) / 2); + int mgrids_num = biggrid_->get_mgrids_num(); + int atom_pair_idx = 0; + for(int i = 0; i < atoms_num; ++i) + { + // only calculate the upper triangle matrix + for(int j = i; j < atoms_num; ++j) + { + int start_idx = mgrids_num; + int end_idx = -1; + for(int mgrid_idx = 0; mgrid_idx < mgrids_num; ++mgrid_idx) + { + if(is_atom_on_mgrid_[i][mgrid_idx] && is_atom_on_mgrid_[j][mgrid_idx]) + { + start_idx = mgrid_idx; + break; + } + } + for(int mgrid_idx = mgrids_num - 1; mgrid_idx >= 0; --mgrid_idx) + { + if(is_atom_on_mgrid_[i][mgrid_idx] && is_atom_on_mgrid_[j][mgrid_idx]) + { + end_idx = mgrid_idx; + break; + } + } + atom_pair_start_end_idx_[atom_pair_idx].first = start_idx; + atom_pair_start_end_idx_[atom_pair_idx].second = end_idx; + atom_pair_idx++; + } + } +} + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.h b/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.h new file mode 100644 index 0000000000..535e3a2f7e --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/phi_operator.h @@ -0,0 +1,129 @@ +#pragma once + +#include +#include +#include +#include +#include "big_grid.h" + +namespace ModuleGint +{ + +/** + * @brief The class phiOperator is used to perform operations on the wave function matrix phi, dphi, etc. + * + * In fact, the variables and functions of this class could be placed in the BigGrid class, but the lifecycle of the BigGrid class is relatively long. + * We do not want the BigGrid to contain too many member variables, as this could lead to excessive memory usage. + * Therefore, we separate this class out, so it can be destroyed after use. + */ +class PhiOperator +{ + public: + // constructor + PhiOperator()=default; + + // set the big grid that the phiOperator is associated with + void set_bgrid(std::shared_ptr biggrid); + + // getter + int get_rows() const {return rows_;}; + int get_cols() const {return cols_;}; + + // get phi of the big grid + // the dimension of phi is num_mgrids * (\sum_{i=0}^{atoms_->size()} atoms_[i]->nw) + void set_phi(double* phi) const; + + // get phi and the gradient of phi of the big grid + // the dimension of phi and dphi is num_mgrids * (\sum_{i=0}^{atoms_->size()} atoms_[i]->nw) + // if you do not need phi, you can set phi to nullptr. + void set_phi_dphi(double* phi, double* dphi_x, double* dphi_y, double* dphi_z) const; + + // get the hessian of the wave function values of the big grid + // the dimension of ddphi is num_mgrids * (\sum_{i=0}^{atoms_->size()} atoms_[i]->nw) + void set_ddphi( + double* ddphi_xx, double* ddphi_xy, double* ddphi_xz, + double* ddphi_yy, double* ddphi_yz, double* ddphi_zz) const; + + void phi_mul_dm( + const double* phi, + const HContainer& DM, + const bool is_symm, double* phi_dm) const; + + void phi_mul_vldr3( + const double* vl, + const double dr3, + const double* phi, + double* result) const; + + void phi_mul_phi_vldr3( + const double* phi, + const double* phi_vldr3, + HContainer* hr) const; + + void phi_dot_phi_dm( + const double* phi, + const double* phi_dm, + double* rho) const; + + void phi_dot_dphi( + const double* phi, + const double* dphi_x, + const double* dphi_y, + const double* dphi_z, + ModuleBase::matrix *fvl) const; + + void phi_dot_dphi_r( + const double* phi, + const double* dphi_x, + const double* dphi_y, + const double* dphi_z, + ModuleBase::matrix *svl) const; + + private: + void init_atom_pair_start_end_idx_(); + + // get the index of the first and the last meshgrid that both atom a and atom b affect + // Note that atom_pair_start_end_idx_ only stores the cases where a <= b, so this function is needed to retrieve the value + const std::pair& get_atom_pair_start_end_idx_(int a, int b) const + { + int x = std::min(a, b); + int y = std::abs(a - b); + return atom_pair_start_end_idx_[(2 * biggrid_->get_atoms_num() - x + 1) * x / 2 + y]; + }; + + // the row number of the phi matrix + // rows_ = biggrid_->get_mgrids_num() + int rows_; + + // the column number of the phi matrix + // cols_ = biggrid_->get_mgrid_phi_len() + int cols_; + + // the local index of the meshgrids + std::vector meshgrids_local_idx_; + + // the big grid that the phi matrix is associated with + std::shared_ptr biggrid_; + + // the relative coordinates of the atoms and the meshgrids + // atoms_relative_coords_[i][j] is the relative coordinate of the jth meshgrid and the ith atom + std::vector> atoms_relative_coords_; + + // record whether the atom affects the meshgrid + // is_atom_on_mgrid_[i][j] = true if the ith atom affects the jth meshgrid, otherwise false + // FIXME,std::vector> is not a efficient data structure, we can use a 1D array to replace it. + std::vector> is_atom_on_mgrid_; + + // the start index of the phi of each atom + std::vector atoms_startidx_; + + // the length of phi of each atom + // atoms_phi_len_[i] = biggrid_->get_atom(i)->get_nw() + // TODO: remove it + std::vector atoms_phi_len_; + + // This data structure is used to store the index of the first and last meshgrid affected by each atom pair + std::vector> atom_pair_start_end_idx_; +}; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/set_ddphi.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/set_ddphi.cpp new file mode 100644 index 0000000000..dd0b77f18f --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/set_ddphi.cpp @@ -0,0 +1,264 @@ +#include "module_base/array_pool.h" +#include "module_base/timer.h" +#include "module_base/ylm.h" +#include "gint_atom.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +template +void GintAtom::set_ddphi( + const std::vector& coords, const int stride, + T* ddphi_xx, T* ddphi_xy, T* ddphi_xz, + T* ddphi_yy, T* ddphi_yz, T* ddphi_zz) const +{ + ModuleBase::timer::tick("GintAtom", "set_ddphi"); + + const int num_mgrids = coords.size(); + + // orb_ does not have the member variable dr_uniform + const double dr_uniform = orb_->PhiLN(0, 0).dr_uniform; + + // store the pointer to reduce repeated address fetching + std::vector p_psi_uniform(atom_->nw); + std::vector p_dpsi_uniform(atom_->nw); + std::vector p_ddpsi_uniform(atom_->nw); + std::vector phi_nr_uniform(atom_->nw); + for (int iw=0; iw< atom_->nw; ++iw) + { + if ( atom_->iw2_new[iw] ) + { + int l = atom_->iw2l[iw]; + int n = atom_->iw2n[iw]; + p_psi_uniform[iw] = orb_->PhiLN(l, n).psi_uniform.data(); + p_dpsi_uniform[iw] = orb_->PhiLN(l, n).dpsi_uniform.data(); + p_ddpsi_uniform[iw] = orb_->PhiLN(l, n).ddpsi_uniform.data(); + phi_nr_uniform[iw] = orb_->PhiLN(l, n).nr_uniform; + } + } + + std::vector rly(std::pow(atom_->nwl + 1, 2)); + ModuleBase::Array_Pool grly(std::pow(atom_->nwl + 1, 2), 3); + // TODO: A better data structure such as a 3D tensor can be used to store dphi + std::vector>> dphi(atom_->nw, std::vector>(6, std::vector(3))); + Vec3d coord1; + ModuleBase::Array_Pool displ(6, 3); + displ[0][0] = 0.0001; // in x direction + displ[1][0] = -0.0001; + displ[2][1] = 0.0001; // in y direction + displ[3][1] = -0.0001; + displ[4][2] = 0.0001; // in z direction + displ[5][2] = -0.0001; + + for(int im = 0; im < num_mgrids; im++) + { + const Vec3d& coord = coords[im]; + // 1e-9 is to avoid division by zero + const double dist = coord.norm() < 1e-9 ? 1e-9 : coord.norm(); + + if(dist > orb_->getRcut()) + { + // if the distance is larger than the cutoff radius, + // the wave function values are all zeros + ModuleBase::GlobalFunc::ZEROS(ddphi_xx + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(ddphi_xy + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(ddphi_xz + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(ddphi_yy + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(ddphi_yz + im * stride, atom_->nw); + ModuleBase::GlobalFunc::ZEROS(ddphi_zz + im * stride, atom_->nw); + continue; + } + + for(int i = 0; i < 6; i++) + { + coord1[0] = coord[0] + displ[i][0]; + coord1[1] = coord[1] + displ[i][1]; + coord1[2] = coord[2] + displ[i][2]; + + // sphereical harmonics + ModuleBase::Ylm::grad_rl_sph_harm(atom_->nwl, coord1[0], coord1[1], coord1[2], rly.data(), grly.get_ptr_2D()); + + const double dist1 = coord1.norm() < 1e-9 ? 1e-9 : coord1.norm(); + + const double position = dist1 / dr_uniform; + const int ip = static_cast(position); + const double x0 = position - ip; + const double x1 = 1.0 - x0; + const double x2 = 2.0 - x0; + const double x3 = 3.0 - x0; + const double x12 = x1 * x2 / 6; + const double x03 = x0 * x3 / 2; + + double tmp, dtmp; + + for(int iw = 0; iw < atom_->nw; ++iw) + { + if(atom_->iw2_new[iw]) + { + auto psi_uniform = p_psi_uniform[iw]; + auto dpsi_uniform = p_dpsi_uniform[iw]; + + if(ip >= phi_nr_uniform[iw] - 4) + { + tmp = dtmp = 0.0; + } + else + { + // use Polynomia Interpolation method to get the + // wave functions + + tmp = x12 * (psi_uniform[ip] * x3 + psi_uniform[ip + 3] * x0) + + x03 * (psi_uniform[ip + 1] * x2 - psi_uniform[ip + 2] * x1); + + dtmp = x12 * (dpsi_uniform[ip] * x3 + dpsi_uniform[ip + 3] * x0) + + x03 * (dpsi_uniform[ip + 1] * x2 - dpsi_uniform[ip + 2] * x1); + } + } + + // get the 'l' of this localized wave function + const int ll = atom_->iw2l[iw]; + const int idx_lm = atom_->iw2_ylm[iw]; + + const double rl = pow_int(dist1, ll); + + // derivative of wave functions with respect to atom positions. + const double tmpdphi_rly = (dtmp - tmp * ll / dist1) / rl * rly[idx_lm] / dist1; + const double tmprl = tmp / rl; + + dphi[iw][i][0] = tmpdphi_rly * coord1[0] + tmprl * grly[idx_lm][0]; + dphi[iw][i][1] = tmpdphi_rly * coord1[1] + tmprl * grly[idx_lm][1]; + dphi[iw][i][2] = tmpdphi_rly * coord1[2] + tmprl * grly[idx_lm][2]; + } // end iw + } // end i + + for(int iw = 0; iw < atom_->nw; iw++) + { + int idx = im * stride + iw; + ddphi_xx[idx] = (dphi[iw][0][0] - dphi[iw][1][0]) / 0.0002; + ddphi_xy[idx] + = ((dphi[iw][2][0] - dphi[iw][3][0]) + (dphi[iw][0][1] - dphi[iw][1][1])) / 0.0004; + ddphi_xz[idx] + = ((dphi[iw][4][0] - dphi[iw][5][0]) + (dphi[iw][0][2] - dphi[iw][1][2])) / 0.0004; + ddphi_yy[idx] = (dphi[iw][2][1] - dphi[iw][3][1]) / 0.0002; + ddphi_yz[idx] + = ((dphi[iw][4][1] - dphi[iw][5][1]) + (dphi[iw][2][2] - dphi[iw][3][2])) / 0.0004; + ddphi_zz[idx] = (dphi[iw][4][2] - dphi[iw][5][2]) / 0.0002; + } + + // else + // // the analytical method for evaluating 2nd derivatives + // // it is not used currently + // { + // // Add it here, but do not run it. If there is a need to run this code + // // in the future, include it in the previous initialization process. + // for (int iw=0; iw< atom->nw; ++iw) + // { + // if ( atom->iw2_new[iw] ) + // { + // it_ddpsi_uniform[iw] = gt.d2phi_u[it*gt.nwmax + iw].data(); + // } + // } + // // End of code addition section. + + // std::vector> hrly; + // ModuleBase::Ylm::grad_rl_sph_harm(ucell.atoms[it].nwl, dr[0], dr[1], dr[2], rly, grly.data()); + // ModuleBase::Ylm::hes_rl_sph_harm(ucell.atoms[it].nwl, dr[0], dr[1], dr[2], hrly); + // const double position = distance / delta_r; + + // const double iq = static_cast(position); + // const int ip = static_cast(position); + // const double x0 = position - iq; + // const double x1 = 1.0 - x0; + // const double x2 = 2.0 - x0; + // const double x3 = 3.0 - x0; + // const double x12 = x1 * x2 / 6; + // const double x03 = x0 * x3 / 2; + + // double tmp, dtmp, ddtmp; + + // for (int iw = 0; iw < atom->nw; ++iw) + // { + // // this is a new 'l', we need 1D orbital wave + // // function from interpolation method. + // if (atom->iw2_new[iw]) + // { + // auto psi_uniform = it_psi_uniform[iw]; + // auto dpsi_uniform = it_dpsi_uniform[iw]; + // auto ddpsi_uniform = it_ddpsi_uniform[iw]; + + // // if ( iq[id] >= philn.nr_uniform-4) + // if (iq >= it_phi_nr_uniform[iw]-4) + // { + // tmp = dtmp = ddtmp = 0.0; + // } + // else + // { + // // use Polynomia Interpolation method to get the + // // wave functions + + // tmp = x12 * (psi_uniform[ip] * x3 + psi_uniform[ip + 3] * x0) + // + x03 * (psi_uniform[ip + 1] * x2 - psi_uniform[ip + 2] * x1); + + // dtmp = x12 * (dpsi_uniform[ip] * x3 + dpsi_uniform[ip + 3] * x0) + // + x03 * (dpsi_uniform[ip + 1] * x2 - dpsi_uniform[ip + 2] * x1); + + // ddtmp = x12 * (ddpsi_uniform[ip] * x3 + ddpsi_uniform[ip + 3] * x0) + // + x03 * (ddpsi_uniform[ip + 1] * x2 - ddpsi_uniform[ip + 2] * x1); + // } + // } // new l is used. + + // // get the 'l' of this localized wave function + // const int ll = atom->iw2l[iw]; + // const int idx_lm = atom->iw2_ylm[iw]; + + // const double rl = pow_int(distance, ll); + // const double r_lp2 =rl * distance * distance; + + // // d/dr (R_l / r^l) + // const double tmpdphi = (dtmp - tmp * ll / distance) / rl; + // const double term1 = ddtmp / r_lp2; + // const double term2 = (2 * ll + 1) * dtmp / r_lp2 / distance; + // const double term3 = ll * (ll + 2) * tmp / r_lp2 / distance / distance; + // const double term4 = tmpdphi / distance; + // const double term5 = term1 - term2 + term3; + + // // hessian of (R_l / r^l) + // const double term_xx = term4 + dr[0] * dr[0] * term5; + // const double term_xy = dr[0] * dr[1] * term5; + // const double term_xz = dr[0] * dr[2] * term5; + // const double term_yy = term4 + dr[1] * dr[1] * term5; + // const double term_yz = dr[1] * dr[2] * term5; + // const double term_zz = term4 + dr[2] * dr[2] * term5; + + // // d/dr (R_l / r^l) * alpha / r + // const double term_1x = dr[0] * term4; + // const double term_1y = dr[1] * term4; + // const double term_1z = dr[2] * term4; + + // p_ddphi_xx[iw] + // = term_xx * rly[idx_lm] + 2.0 * term_1x * grly[idx_lm][0] + tmp / rl * hrly[idx_lm][0]; + // p_ddphi_xy[iw] = term_xy * rly[idx_lm] + term_1x * grly[idx_lm][1] + term_1y * grly[idx_lm][0] + // + tmp / rl * hrly[idx_lm][1]; + // p_ddphi_xz[iw] = term_xz * rly[idx_lm] + term_1x * grly[idx_lm][2] + term_1z * grly[idx_lm][0] + // + tmp / rl * hrly[idx_lm][2]; + // p_ddphi_yy[iw] + // = term_yy * rly[idx_lm] + 2.0 * term_1y * grly[idx_lm][1] + tmp / rl * hrly[idx_lm][3]; + // p_ddphi_yz[iw] = term_yz * rly[idx_lm] + term_1y * grly[idx_lm][2] + term_1z * grly[idx_lm][1] + // + tmp / rl * hrly[idx_lm][4]; + // p_ddphi_zz[iw] + // = term_zz * rly[idx_lm] + 2.0 * term_1z * grly[idx_lm][2] + tmp / rl * hrly[idx_lm][5]; + + // } // iw + // } // end if + } + + ModuleBase::timer::tick("GintAtom", "set_ddphi"); +} + +// explicit instantiation +template void GintAtom::set_ddphi(const std::vector& coords, const int stride, + double* ddphi_xx, double* ddphi_xy, double* ddphi_xz, + double* ddphi_yy, double* ddphi_yz, double* ddphi_zz) const; + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.cpp b/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.cpp new file mode 100644 index 0000000000..7818b3b08a --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.cpp @@ -0,0 +1,29 @@ +#include "unitcell_info.h" +#include "gint_helper.h" + +namespace ModuleGint +{ + +UnitCellInfo::UnitCellInfo( + Vec3d unitcell_vec1, + Vec3d unitcell_vec2, + Vec3d unitcell_vec3, + int nbx, int nby, int nbz, + int nmx, int nmy, int nmz) + : unitcell_vec1_(unitcell_vec1), + unitcell_vec2_(unitcell_vec2), + unitcell_vec3_(unitcell_vec3), + nbx_(nbx), nby_(nby), nbz_(nbz), nbxyz_(nbx*nby*nbz), + nmx_(nmx), nmy_(nmy), nmz_(nmz), nmxyz_(nmx*nmy*nmz) + { + // initialize the biggrid_info_ + biggrid_info_ = std::make_shared( + unitcell_vec1_ / static_cast(nbx), + unitcell_vec2_ / static_cast(nby), + unitcell_vec3_ / static_cast(nbz), + nmx/nbx, nmy/nby, nmz/nbz); + + meshgrid_info_ = biggrid_info_->get_mgrid_info(); + } + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.h b/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.h new file mode 100644 index 0000000000..868f51d830 --- /dev/null +++ b/source/module_hamilt_lcao/module_gint/temp_gint/unitcell_info.h @@ -0,0 +1,172 @@ +#pragma once + +#include +#include +#include "biggrid_info.h" +#include "gint_helper.h" +#include "gint_type.h" + +namespace ModuleGint +{ + +class UnitCellInfo +{ + public: + // constructor + UnitCellInfo( + Vec3d unitcell_vec1, + Vec3d unitcell_vec2, + Vec3d unitcell_vec3, + int nbx, int nby, int nbz, + int nmx, int nmy, int nmz); + + // getter functions + int get_nbx() const { return nbx_; }; + int get_nby() const { return nby_; }; + int get_nbz() const { return nbz_; }; + int get_bgrids_num() const { return nbxyz_; }; + int get_nmx() const { return nmx_; }; + int get_nmy() const { return nmy_; }; + int get_nmz() const { return nmz_; }; + int get_mgrids_num() const { return nmxyz_; }; + std::shared_ptr get_bgrid_info() const { return biggrid_info_; }; + std::shared_ptr get_mgrid_info() const { return meshgrid_info_; }; + + //==================================================================== + // functions related to the big grid + //==================================================================== + + // transform the 1D index of a big grid in the unit cell to the 3D index + Vec3i bgrid_idx_1Dto3D(const int index_1d) const + { + return index1Dto3D(index_1d, nbx_, nby_, nbz_); + }; + + // transform the 3D index of a biggrid in the unit cell to the 1D index + int bgrid_idx_3Dto1D(const Vec3i index_3d) const + { + return index3Dto1D(index_3d.x, index_3d.y, index_3d.z, nbx_, nby_, nbz_); + }; + + // get the cartesian coordinate of a big grid in the unit cell from the 3D index + Vec3d get_bgrid_coord(Vec3i index_3d) const + { + return biggrid_info_->get_cartesian_coord(index_3d); + }; + + // get the cartesian coordinate of a big grid in the unit cell from the 1D index + Vec3d get_bgrid_coord(int index_1d) const + { + return get_bgrid_coord(bgrid_idx_1Dto3D(index_1d)); + }; + + // get the 3D index of a big grid in the unit cell from the cartesian coordinate + Vec3i get_bgrid_idx_3d(const Vec3d coord) const + { + Vec3d direct_coord = biggrid_info_->get_direct_coord(coord); + return Vec3i( + static_cast(floor(direct_coord.x)), + static_cast(floor(direct_coord.y)), + static_cast(floor(direct_coord.z))); + }; + + // Get the relative Cartesian coordinates of big grid A relative to big grid B + // returned vector = coordinates of point A - coordinates of point B + // this function is more efficient than obtaining two 3D coordinates separately + // through two 3D indices and then subtracting them + Vec3d get_relative_coord(Vec3i index_3d_a, Vec3i index_3d_b) const + { + return get_bgrid_coord(index_3d_a - index_3d_b); + }; + + // get the extended unitcell index of a big grid + Vec3i get_unitcell_idx(const Vec3i index_3d) const + { + return Vec3i(floor_div(index_3d.x, nbx_), + floor_div(index_3d.y, nby_), + floor_div(index_3d.z, nbz_)); + }; + + // map the extended big grid index to the big grid index in unitcell + Vec3i map_ext_idx_to_ucell(const Vec3i index_3d) const + { + return Vec3i(index_3d.x - floor_div(index_3d.x, nbx_) * nbx_, + index_3d.y - floor_div(index_3d.y, nby_) * nby_, + index_3d.z - floor_div(index_3d.z, nbz_) * nbz_); + }; + + + //==================================================================== + // functions related to the meshgrid + //==================================================================== + + // transform the 1D index of a meshgrid in the unit cell to the 3D index + Vec3i mgrid_idx_1Dto3D(const int index_1d) const + { + return index1Dto3D(index_1d, nmx_, nmy_, nmz_); + } + + // transform the 3D index of a meshgrid in the unit cell to the 1D index + int mgrid_idx_3Dto1D(const Vec3i index_3d) const + { + return index3Dto1D(index_3d.x, index_3d.y, index_3d.z, nmx_, nmy_, nmz_); + } + + // get the cartesian coordinate of a meshgrid in the unit cell from the 3D index + Vec3d get_mgrid_coord(Vec3i index_3d) const + { + return meshgrid_info_->get_cartesian_coord(index_3d); + }; + + // get the cartesian coordinate of a meshgrid in the unit cell from the 1D index + Vec3d get_mgrid_coord(int index_1d) const + { + return get_mgrid_coord(mgrid_idx_1Dto3D(index_1d)); + } + + private: + // basis vectors of the unit cell + Vec3d unitcell_vec1_; + Vec3d unitcell_vec2_; + Vec3d unitcell_vec3_; + + //==================================================================== + // member variables related to the Big Grid + //==================================================================== + + // the number of big cells along the first lattice vector + int nbx_; + + // the number of big cells along the second lattice vector + int nby_; + + // the number of big cells along the third lattice vector + int nbz_; + + // the total number of big cells + int nbxyz_; + + // basic attributes of the big grid + std::shared_ptr biggrid_info_; + + std::shared_ptr meshgrid_info_; + + //==================================================================== + // member variables related to meshgrid + //==================================================================== + + // the number of meshgrids along the first lattice vector + int nmx_; + + // the number of meshgrids along the second lattice vector + int nmy_; + + // the number of meshgrids along the third lattice vector + int nmz_; + + // the total number of meshgrids in the unitcell + int nmxyz_; + +}; + +} // namespace ModuleGint \ No newline at end of file diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp index dc0398e8c9..abe56d71c3 100644 --- a/source/module_rdmft/update_state_rdmft.cpp +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -8,7 +8,7 @@ #include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_elecstate/module_charge/symmetry_rho.h" -#include "module_elecstate/elecstate_lcao_cal_tau.h" + namespace rdmft { @@ -118,7 +118,7 @@ void RDMFT::update_charge(UnitCell& ucell) // } // Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau); // GG->cal_gint(&inout1); - elecstate::lcao_cal_tau_gamma(GG, charge); + this->pelec->cal_tau(wfc); } charge->renormalize_rho(); @@ -148,7 +148,7 @@ void RDMFT::update_charge(UnitCell& ucell) // } // Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau); // GK->cal_gint(&inout1); - elecstate::lcao_cal_tau_k(GK, charge); + this->pelec->cal_tau(wfc); } charge->renormalize_rho();