diff --git a/source/module_basis/module_ao/ORB_nonlocal_lm.h b/source/module_basis/module_ao/ORB_nonlocal_lm.h index 31ed209f14..b3504065ae 100644 --- a/source/module_basis/module_ao/ORB_nonlocal_lm.h +++ b/source/module_basis/module_ao/ORB_nonlocal_lm.h @@ -45,6 +45,9 @@ class Numerical_Nonlocal_Lm const double& getKpoint(const int &ik) const { return this->k_radial[ik]; } const double* getBeta_k() const { return this->beta_k; } const double& getBeta_k(const int &ik) const { return this->beta_k[ik]; } + + const int& getNk() const { return nk; } + const double& getDruniform() const { return dr_uniform; } // enables deep copy Numerical_Nonlocal_Lm& operator= (const Numerical_Nonlocal_Lm& nol ); diff --git a/source/module_io/cal_r_overlap_R.cpp b/source/module_io/cal_r_overlap_R.cpp index dec4cbabe8..384b61677c 100644 --- a/source/module_io/cal_r_overlap_R.cpp +++ b/source/module_io/cal_r_overlap_R.cpp @@ -5,6 +5,7 @@ #include "module_base/timer.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/mathzone_add1.h" cal_r_overlap_R::cal_r_overlap_R() { @@ -224,6 +225,230 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const LCAO_Orbitals& orb) } } +void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const LCAO_Orbitals& orb) +{ + const InfoNonlocal& infoNL_ = GlobalC::ucell.infoNL; + + int orb_r_ntype = 0; + int mat_Nr = orb.Phi[0].PhiLN(0, 0).getNr(); + int count_Nr = 0; + + orbs.resize(orb.get_ntype()); + for (int T = 0; T < orb.get_ntype(); ++T) + { + count_Nr = orb.Phi[T].PhiLN(0, 0).getNr(); + if (count_Nr > mat_Nr) + { + mat_Nr = count_Nr; + orb_r_ntype = T; + } + + orbs[T].resize(orb.Phi[T].getLmax() + 1); + for (int L = 0; L <= orb.Phi[T].getLmax(); ++L) + { + orbs[T][L].resize(orb.Phi[T].getNchi(L)); + for (int N = 0; N < orb.Phi[T].getNchi(L); ++N) + { + const auto& orb_origin = orb.Phi[T].PhiLN(L, N); + orbs[T][L][N].set_orbital_info(orb_origin.getLabel(), + orb_origin.getType(), + orb_origin.getL(), + orb_origin.getChi(), + orb_origin.getNr(), + orb_origin.getRab(), + orb_origin.getRadial(), + Numerical_Orbital_Lm::Psi_Type::Psi, + orb_origin.getPsi(), + static_cast(orb_origin.getNk() * kmesh_times) | 1, + orb_origin.getDk(), + orb_origin.getDruniform(), + false, + true, + PARAM.inp.cal_force); + } + } + } + + orb_r.set_orbital_info(orbs[orb_r_ntype][0][0].getLabel(), // atom label + orb_r_ntype, // atom type + 1, // angular momentum L + 1, // number of orbitals of this L , just N + orbs[orb_r_ntype][0][0].getNr(), // number of radial mesh + orbs[orb_r_ntype][0][0].getRab(), // the mesh interval in radial mesh + orbs[orb_r_ntype][0][0].getRadial(), // radial mesh value(a.u.) + Numerical_Orbital_Lm::Psi_Type::Psi, + orbs[orb_r_ntype][0][0].getRadial(), // radial wave function + orbs[orb_r_ntype][0][0].getNk(), + orbs[orb_r_ntype][0][0].getDk(), + orbs[orb_r_ntype][0][0].getDruniform(), + false, + true, + PARAM.inp.cal_force); + + orbs_nonlocal.resize(orb.get_ntype()); + for (int T = 0; T < orb.get_ntype(); ++T) + { + const int nproj = infoNL_.nproj[T]; + orbs_nonlocal[T].resize(nproj); + for (int ip = 0; ip < nproj; ip++) + { + int nr = infoNL_.Beta[T].Proj[ip].getNr(); + double dr_uniform = 0.01; + int nr_uniform = static_cast((infoNL_.Beta[T].Proj[ip].getRadial(nr-1) - infoNL_.Beta[T].Proj[ip].getRadial(0))/dr_uniform) + 1; + double* rad = new double[nr_uniform]; + double* rab = new double[nr_uniform]; + for (int ir = 0; ir < nr_uniform; ir++) + { + rad[ir] = ir*dr_uniform; + rab[ir] = dr_uniform; + } + double* y2 = new double[nr]; + double* Beta_r_uniform = new double[nr_uniform]; + double* dbeta_uniform = new double[nr_uniform]; + ModuleBase::Mathzone_Add1::SplineD2(infoNL_.Beta[T].Proj[ip].getRadial(), infoNL_.Beta[T].Proj[ip].getBeta_r(), nr, 0.0, 0.0, y2); + ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation( + infoNL_.Beta[T].Proj[ip].getRadial(), + infoNL_.Beta[T].Proj[ip].getBeta_r(), + y2, + nr, + rad, + nr_uniform, + Beta_r_uniform, + dbeta_uniform + ); + + // linear extrapolation at the zero point + if (infoNL_.Beta[T].Proj[ip].getRadial(0) > 1e-10) + { + double slope = (infoNL_.Beta[T].Proj[ip].getBeta_r(1) - infoNL_.Beta[T].Proj[ip].getBeta_r(0)) / (infoNL_.Beta[T].Proj[ip].getRadial(1) - infoNL_.Beta[T].Proj[ip].getRadial(0)); + Beta_r_uniform[0] = infoNL_.Beta[T].Proj[ip].getBeta_r(0) - slope * infoNL_.Beta[T].Proj[ip].getRadial(0); + } + + orbs_nonlocal[T][ip].set_orbital_info(infoNL_.Beta[T].getLabel(), + infoNL_.Beta[T].getType(), + infoNL_.Beta[T].Proj[ip].getL(), + 1, + nr_uniform, + rab, + rad, + Numerical_Orbital_Lm::Psi_Type::Psi, + Beta_r_uniform, + static_cast(infoNL_.Beta[T].Proj[ip].getNk() * kmesh_times) | 1, + infoNL_.Beta[T].Proj[ip].getDk(), + infoNL_.Beta[T].Proj[ip].getDruniform(), + false, + true, + PARAM.inp.cal_force); + + delete [] rad; + delete [] rab; + delete [] y2; + delete [] Beta_r_uniform; + delete [] dbeta_uniform; + } + } + + for (int TA = 0; TA < orb.get_ntype(); ++TA) + { + for (int TB = 0; TB < orb.get_ntype(); ++TB) + { + for (int LA = 0; LA <= orb.Phi[TA].getLmax(); ++LA) + { + for (int NA = 0; NA < orb.Phi[TA].getNchi(LA); ++NA) + { + for (int ip = 0; ip < infoNL_.nproj[TB]; ip++) + { + center2_orb11_nonlocal[TA][TB][LA][NA].insert( + std::make_pair(ip, Center2_Orb::Orb11(orbs[TA][LA][NA], orbs_nonlocal[TB][ip], psb_, MGT))); + } + } + } + } + } + + for (int TA = 0; TA < orb.get_ntype(); ++TA) + { + for (int TB = 0; TB < orb.get_ntype(); ++TB) + { + for (int LA = 0; LA <= orb.Phi[TA].getLmax(); ++LA) + { + for (int NA = 0; NA < orb.Phi[TA].getNchi(LA); ++NA) + { + for (int ip = 0; ip < infoNL_.nproj[TB]; ip++) + { + center2_orb21_r_nonlocal[TA][TB][LA][NA].insert(std::make_pair( + ip, + Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs_nonlocal[TB][ip], psb_, MGT))); + } + } + } + } + } + + for (auto& co1: center2_orb11_nonlocal) + { + for (auto& co2: co1.second) + { + for (auto& co3: co2.second) + { + for (auto& co4: co3.second) + { + for (auto& co5: co4.second) + { + co5.second.init_radial_table(); + } + } + } + } + } + + for (auto& co1: center2_orb21_r_nonlocal) + { + for (auto& co2: co1.second) + { + for (auto& co3: co2.second) + { + for (auto& co4: co3.second) + { + for (auto& co5: co4.second) + { + co5.second.init_radial_table(); + } + } + } + } + } + + iw2it.resize(PARAM.globalv.nlocal); + iw2ia.resize(PARAM.globalv.nlocal); + iw2iL.resize(PARAM.globalv.nlocal); + iw2iN.resize(PARAM.globalv.nlocal); + iw2im.resize(PARAM.globalv.nlocal); + + int iw = 0; + for (int it = 0; it < GlobalC::ucell.ntype; it++) + { + for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++) + { + for (int iL = 0; iL < GlobalC::ucell.atoms[it].nwl + 1; iL++) + { + for (int iN = 0; iN < GlobalC::ucell.atoms[it].l_nchi[iL]; iN++) + { + for (int im = 0; im < (2 * iL + 1); im++) + { + iw2it[iw] = it; + iw2ia[iw] = ia; + iw2iL[iw] = iL; + iw2iN[iw] = iN; + iw2im[iw] = im; + iw++; + } + } + } + } + } +} + void cal_r_overlap_R::init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) { ModuleBase::TITLE("cal_r_overlap_R", "init"); @@ -236,6 +461,21 @@ void cal_r_overlap_R::init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb ModuleBase::timer::tick("cal_r_overlap_R", "init"); return; } + +void cal_r_overlap_R::init_nonlocal(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) +{ + ModuleBase::TITLE("cal_r_overlap_R", "init_nonlocal"); + ModuleBase::timer::tick("cal_r_overlap_R", "init_nonlocal"); + this->ParaV = &pv; + + initialize_orb_table(orb); + construct_orbs_and_nonlocal_and_orb_r(orb); + + ModuleBase::timer::tick("cal_r_overlap_R", "init_nonlocal"); + return; +} + + ModuleBase::Vector3 cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vector3& R1, const int& T1, const int& L1, @@ -279,6 +519,72 @@ ModuleBase::Vector3 cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vec return temp_prp; } + +void cal_r_overlap_R::get_psi_r_beta(std::vector>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2) +{ + ModuleBase::Vector3 origin_point(0.0, 0.0, 0.0); + double factor = sqrt(ModuleBase::FOUR_PI / 3.0); + const ModuleBase::Vector3& distance = R2 - R1; + const InfoNonlocal& infoNL_ = GlobalC::ucell.infoNL; + const int nproj = infoNL_.nproj[T2]; + if (nproj == 0) + { + nlm.resize(1); + return; + } + + int natomwfc = 0; + for (int ip = 0; ip < nproj; ip++) + { + const int L2 = infoNL_.Beta[T2].Proj[ip].getL(); // mohan add 2021-05-07 + natomwfc += 2 * L2 + 1; + } + nlm.resize(natomwfc); + + int index = 0; + for (int ip = 0; ip < nproj; ip++) + { + const int L2 = infoNL_.Beta[T2].Proj[ip].getL(); + for (int m2 = 0; m2 < 2 * L2 + 1; m2++) + { + double overlap_o + = center2_orb11_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, distance, m1, m2); + + double overlap_x = -1 * factor + * center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, + distance, + m1, + 1, + m2); // m = 1 + + double overlap_y = -1 * factor + * center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, + distance, + m1, + 2, + m2); // m = -1 + + double overlap_z = factor + * center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, + distance, + m1, + 0, + m2); // m = 0 + + nlm[index] = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + R1 * overlap_o; + index++; + } + } +} + + void cal_r_overlap_R::out_rR(const int& istep) { ModuleBase::TITLE("cal_r_overlap_R", "out_rR"); diff --git a/source/module_io/cal_r_overlap_R.h b/source/module_io/cal_r_overlap_R.h index 27006123c5..e9e02a14f3 100644 --- a/source/module_io/cal_r_overlap_R.h +++ b/source/module_io/cal_r_overlap_R.h @@ -31,7 +31,8 @@ class cal_r_overlap_R bool binary = false; void init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); - ModuleBase::Vector3 get_psi_r_psi( + void init_nonlocal(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); + ModuleBase::Vector3 get_psi_r_psi( const ModuleBase::Vector3& R1, const int& T1, const int& L1, @@ -43,6 +44,16 @@ class cal_r_overlap_R const int& m2, const int& N2 ); + void get_psi_r_beta( + std::vector>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2 + ); void out_rR(const int& istep); void out_rR_other(const int& istep, const std::set>& output_R_coor); @@ -50,6 +61,7 @@ class cal_r_overlap_R private: void initialize_orb_table(const LCAO_Orbitals& orb); void construct_orbs_and_orb_r(const LCAO_Orbitals& orb); + void construct_orbs_and_nonlocal_and_orb_r(const LCAO_Orbitals& orb); std::vector iw2ia; std::vector iw2iL;