From 48cb0c1285685d03c2133d93304909b9c92084c6 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 11:10:23 +0800 Subject: [PATCH 01/10] update ucell in force_gamma/force_k --- .../module_hamilt_lcao/hamilt_lcaodft/FORCE.h | 3 ++- .../hamilt_lcaodft/FORCE_gamma.cpp | 9 +++---- .../hamilt_lcaodft/FORCE_k.cpp | 24 ++++++++++--------- 3 files changed, 20 insertions(+), 16 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h index ecd6d11b6a..171c94a5d8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h @@ -72,7 +72,8 @@ class Force_LCAO Record_adj* ra = nullptr); // get the ds, dt, dvnl. - void allocate(const Parallel_Orbitals& pv, + void allocate(const UnitCell& ucell, + const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index dcf34c95b1..c7aad83123 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -16,7 +16,8 @@ #include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h" template <> -void Force_LCAO::allocate(const Parallel_Orbitals& pv, +void Force_LCAO::allocate(const UnitCell& ucell, + const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, @@ -75,7 +76,7 @@ void Force_LCAO::allocate(const Parallel_Orbitals& pv, 'S', cal_deri, PARAM.inp.cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -99,7 +100,7 @@ void Force_LCAO::allocate(const Parallel_Orbitals& pv, 'T', cal_deri, PARAM.inp.cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -205,7 +206,7 @@ void Force_LCAO::ftable(const bool isforce, // allocate DSloc_x, DSloc_y, DSloc_z // allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z - this->allocate(pv, fsr, two_center_bundle, orb); + this->allocate(ucell,pv, fsr, two_center_bundle, orb); const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z }; const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 }; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index 7f6d4d1c1c..f00b089a75 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -26,7 +26,8 @@ #endif template <> -void Force_LCAO>::allocate(const Parallel_Orbitals& pv, +void Force_LCAO>::allocate(const UnitCell& ucell, + const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan add 2024-06-15 const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, @@ -93,7 +94,7 @@ void Force_LCAO>::allocate(const Parallel_Orbitals& pv, 'S', cal_deri, PARAM.inp.cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -124,7 +125,7 @@ void Force_LCAO>::allocate(const Parallel_Orbitals& pv, 'T', cal_deri, PARAM.inp.cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -143,7 +144,7 @@ void Force_LCAO>::allocate(const Parallel_Orbitals& pv, 'S', cal_deri, PARAM.inp.cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -210,19 +211,19 @@ void Force_LCAO>::finish_ftable(ForceStressArrays& fsr) // test = new double[PARAM.globalv.nlocal * PARAM.globalv.nlocal]; // ModuleBase::GlobalFunc::ZEROS(test, PARAM.globalv.nlocal * PARAM.globalv.nlocal); // -// for (int T1 = 0; T1 < GlobalC::ucell.ntype; T1++) +// for (int T1 = 0; T1 < ucell.ntype; T1++) // { -// Atom* atom1 = &GlobalC::ucell.atoms[T1]; +// Atom* atom1 = &ucell.atoms[T1]; // for (int I1 = 0; I1 < atom1->na; I1++) // { -// // const int iat = GlobalC::ucell.itia2iat(T1,I1); -// const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); +// // const int iat = ucell.itia2iat(T1,I1); +// const int start1 = ucell.itiaiw2iwt(T1, I1, 0); // for (int cb = 0; cb < RA.na_each[ca]; cb++) // { // const int T2 = RA.info[ca][cb][3]; // const int I2 = RA.info[ca][cb][4]; -// Atom* atom2 = &GlobalC::ucell.atoms[T2]; -// const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); +// Atom* atom2 = &ucell.atoms[T2]; +// const int start2 = ucell.itiaiw2iwt(T2, I2, 0); // // for (int jj = 0; jj < atom1->nw; jj++) // { @@ -296,7 +297,8 @@ void Force_LCAO>::ftable(const bool isforce, elecstate::DensityMatrix, double>* dm = dynamic_cast>*>(pelec)->get_DM(); - this->allocate(pv, + this->allocate(ucell, + pv, fsr, // mohan add 2024-06-16 two_center_bundle, orb, From 9e2be34bc9173ef3e2ccfd90bbd1d28e553416e0 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 11:37:41 +0800 Subject: [PATCH 02/10] change GlobalC::ucell in the force_stress --- source/module_esolver/esolver_ks_lcao.cpp | 1 + .../hamilt_lcaodft/FORCE_STRESS.cpp | 156 ++++++++++-------- .../hamilt_lcaodft/FORCE_STRESS.h | 12 +- 3 files changed, 93 insertions(+), 76 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index a1bdd0cde2..9331c206ff 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -294,6 +294,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for PARAM.inp.cal_stress, PARAM.inp.test_force, PARAM.inp.test_stress, + ucell, this->pv, this->pelec, this->psi, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 7423ade290..cdb0e13290 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -35,6 +35,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const bool isstress, const bool istestf, const bool istests, + const UnitCell& ucell, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -63,7 +64,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, return; } - const int nat = GlobalC::ucell.nat; + const int nat = ucell.nat; ForceStressArrays fsr; // mohan add 2024-06-15 @@ -93,7 +94,8 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, fscc.create(nat, 3); // calculate basic terms in Force, same method with PW base - this->calForcePwPart(fvl_dvl, + this->calForcePwPart(ucell, + fvl_dvl, fewalds, fcc, fscc, @@ -138,7 +140,8 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, svnl_dalpha.create(3, 3); #endif // calculate basic terms in Stress, similar method with PW base - this->calStressPwPart(sigmadvl, + this->calStressPwPart(ucell, + sigmadvl, sigmahar, sigmaewa, sigmacc, @@ -154,6 +157,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, this->integral_part(PARAM.globalv.gamma_only_local, isforce, isstress, + ucell, fsr, pelec, psi, @@ -181,7 +185,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, nullptr, kv.kvec_d, nullptr, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get() @@ -205,7 +209,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, nullptr, kv.kvec_d, nullptr, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get() @@ -228,14 +232,14 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // jiyy add 2019-05-18, update 2021-05-02 ModuleBase::matrix force_vdw; ModuleBase::matrix stress_vdw; - auto vdw_solver = vdw::make_vdw(GlobalC::ucell, PARAM.inp); + auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); if (vdw_solver != nullptr) { if (isforce) { force_vdw.create(nat, 3); const std::vector>& force_vdw_temp = vdw_solver->get_force(); - for (int iat = 0; iat < GlobalC::ucell.nat; ++iat) + for (int iat = 0; iat < ucell.nat; ++iat) { force_vdw(iat, 0) = force_vdw_temp[iat].x; force_vdw(iat, 1) = force_vdw_temp[iat].y; @@ -253,7 +257,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.efield_flag && isforce) { fefield.create(nat, 3); - elecstate::Efield::compute_force(GlobalC::ucell, fefield); + elecstate::Efield::compute_force(ucell, fefield); } //! atomic forces from E-field of rt-TDDFT @@ -261,7 +265,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.esolver_type == "TDDFT" && isforce) { fefield_tddft.create(nat, 3); - elecstate::Efield::compute_force(GlobalC::ucell, fefield_tddft); + elecstate::Efield::compute_force(ucell, fefield_tddft); } //! atomic forces from gate field @@ -269,7 +273,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.gate_flag && isforce) { fgate.create(nat, 3); - elecstate::Gatefield::compute_force(GlobalC::ucell, fgate); + elecstate::Gatefield::compute_force(ucell, fgate); } //! atomic forces from implicit solvation model @@ -277,7 +281,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (PARAM.inp.imp_sol && isforce) { fsol.create(nat, 3); - GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, rhopw, nlpp.vloc, fsol); + GlobalC::solvent_model.cal_force_sol(ucell, rhopw, nlpp.vloc, fsol); } //! atomic forces from DFT+U (Quxin version) @@ -308,7 +312,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, hamilt::DFTU> tmp_dftu(nullptr, // HK and SK are not used for force&stress kv.kvec_d, nullptr, // HR are not used for force&stress - GlobalC::ucell, + ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), @@ -336,7 +340,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, nullptr, kv.kvec_d, nullptr, - GlobalC::ucell, + ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs() @@ -490,7 +494,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // pengfei 2016-12-20 if (ModuleSymmetry::Symmetry::symm_flag == 1) { - this->forceSymmetry(fcs, symm); + this->forceSymmetry(ucell,fcs, symm); } #ifdef __DEEPKS @@ -500,7 +504,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const std::string file_ftot = PARAM.globalv.global_out_dir + "deepks_ftot.npy"; LCAO_deepks_io::save_npy_f(fcs, file_ftot, - GlobalC::ucell.nat, + ucell.nat, GlobalV::MY_RANK); // Ty/Bohr, F_tot if (PARAM.inp.deepks_scf) @@ -508,7 +512,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; LCAO_deepks_io::save_npy_f(fcs - GlobalC::ld.F_delta, file_fbase, - GlobalC::ucell.nat, + ucell.nat, GlobalV::MY_RANK); // Ry/Bohr, F_base if (!PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now @@ -517,7 +521,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, { const std::vector>& dm_gamma = dynamic_cast*>(pelec)->get_DM()->get_DMK_vector(); - GlobalC::ld.cal_gdmx(dm_gamma[0], GlobalC::ucell, orb, GlobalC::GridD, isstress); + GlobalC::ld.cal_gdmx(dm_gamma[0], ucell, orb, GlobalC::GridD, isstress); } else { @@ -527,7 +531,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, ->get_DMK_vector(); GlobalC::ld.cal_gdmx_k(dm_k, - GlobalC::ucell, + ucell, orb, GlobalC::GridD, kv.get_nks(), @@ -536,16 +540,16 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_gdmx(GlobalC::ucell.nat); + GlobalC::ld.check_gdmx(ucell.nat); } - GlobalC::ld.cal_gvx(GlobalC::ucell.nat); + GlobalC::ld.cal_gvx(ucell.nat); if (PARAM.inp.deepks_out_unittest) { - GlobalC::ld.check_gvx(GlobalC::ucell.nat); + GlobalC::ld.check_gvx(ucell.nat); } - LCAO_deepks_io::save_npy_gvx(GlobalC::ucell.nat, + LCAO_deepks_io::save_npy_gvx(ucell.nat, GlobalC::ld.des_per_atom, GlobalC::ld.gvx_tensor, PARAM.globalv.global_out_dir, @@ -557,7 +561,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; LCAO_deepks_io::save_npy_f(fcs, file_fbase, - GlobalC::ucell.nat, + ucell.nat, GlobalV::MY_RANK); // no scf, F_base=F_tot } } @@ -587,20 +591,20 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // regular force terms test. //----------------------------- // this->print_force("OVERLAP FORCE",foverlap,1,ry); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "OVERLAP FORCE", foverlap, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "OVERLAP FORCE", foverlap, false); // this->print_force("TVNL_DPHI force",ftvnl_dphi,PARAM.inp.test_force); // this->print_force("VNL_DBETA force",fvnl_dbeta,PARAM.inp.test_force); // this->print_force("T_VNL FORCE",ftvnl,1,ry); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "T_VNL FORCE", ftvnl, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "VL_dPHI FORCE", fvl_dphi, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "T_VNL FORCE", ftvnl, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "VL_dPHI FORCE", fvl_dphi, false); // this->print_force("VL_dPHI FORCE",fvl_dphi,1,ry); // this->print_force("VL_dVL FORCE",fvl_dvl,1,ry); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "VL_dVL FORCE", fvl_dvl, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EWALD FORCE", fewalds, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "VL_dVL FORCE", fvl_dvl, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "EWALD FORCE", fewalds, false); // this->print_force("VLOCAL FORCE",fvlocal,PARAM.inp.test_force); // this->print_force("EWALD FORCE",fewalds,1,ry); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "NLCC FORCE", fcc, false); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "SCC FORCE", fscc, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "NLCC FORCE", fcc, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "SCC FORCE", fscc, false); // this->print_force("NLCC FORCE",fcc,1,ry); // this->print_force("SCC FORCE",fscc,1,ry); //------------------------------- @@ -608,42 +612,42 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, //------------------------------- if (PARAM.inp.efield_flag) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD FORCE", fefield, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD FORCE", fefield, false); // this->print_force("EFIELD FORCE",fefield,1,ry); } if (PARAM.inp.esolver_type == "TDDFT") { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "EFIELD_TDDFT FORCE", fefield_tddft, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD_TDDFT FORCE", fefield_tddft, false); // this->print_force("EFIELD_TDDFT FORCE",fefield_tddft,1,ry); } if (PARAM.inp.gate_flag) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "GATEFIELD FORCE", fgate, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "GATEFIELD FORCE", fgate, false); // this->print_force("GATEFIELD FORCE",fgate,1,ry); } if (PARAM.inp.imp_sol) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "IMP_SOL FORCE", fsol, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE", fsol, false); // this->print_force("IMP_SOL FORCE",fsol,1,ry); } if (vdw_solver != nullptr) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "VDW FORCE", force_vdw, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "VDW FORCE", force_vdw, false); // this->print_force("VDW FORCE",force_vdw,1,ry); } if (PARAM.inp.dft_plus_u) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "DFT+U FORCE", force_dftu, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "DFT+U FORCE", force_dftu, false); } if (PARAM.inp.sc_mag_switch) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "DeltaSpin FORCE", force_dspin, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "DeltaSpin FORCE", force_dspin, false); } #ifdef __DEEPKS // caoyu add 2021-06-03 if (PARAM.inp.deepks_scf) { - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "DeePKS FORCE", GlobalC::ld.F_delta, true); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "DeePKS FORCE", GlobalC::ld.F_delta, true); // this->print_force("DeePKS FORCE", GlobalC::ld.F_delta, 1, ry); } #endif @@ -652,13 +656,13 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, GlobalV::ofs_running << std::setiosflags(std::ios::left); // this->printforce_total(ry, istestf, fcs); - ModuleIO::print_force(GlobalV::ofs_running, GlobalC::ucell, "TOTAL-FORCE (eV/Angstrom)", fcs, false); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "TOTAL-FORCE (eV/Angstrom)", fcs, false); if (istestf) { GlobalV::ofs_running << "\n FORCE INVALID TABLE." << std::endl; GlobalV::ofs_running << " " << std::setw(8) << "atom" << std::setw(5) << "x" << std::setw(5) << "y" << std::setw(5) << "z" << std::endl; - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int iat = 0; iat < ucell.nat; iat++) { GlobalV::ofs_running << " " << std::setw(8) << iat; for (int i = 0; i < 3; i++) @@ -719,7 +723,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } if (ModuleSymmetry::Symmetry::symm_flag == 1) { - symm->symmetrize_mat3(scs, GlobalC::ucell.lat); + symm->symmetrize_mat3(scs, ucell.lat); } // end symmetry #ifdef __DEEPKS @@ -728,14 +732,14 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, const std::string file_s = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; LCAO_deepks_io::save_npy_s(scs, file_s, - GlobalC::ucell.omega, + ucell.omega, GlobalV::MY_RANK); // change to energy unit Ry when printing, S_base; } if (PARAM.inp.deepks_scf) { if (ModuleSymmetry::Symmetry::symm_flag == 1) { - symm->symmetrize_mat3(svnl_dalpha, GlobalC::ucell.lat); + symm->symmetrize_mat3(svnl_dalpha, ucell.lat); } // end symmetry for (int i = 0; i < 3; i++) { @@ -751,7 +755,7 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, LCAO_deepks_io::save_npy_s( scs, file_s, - GlobalC::ucell.omega, + ucell.omega, GlobalV::MY_RANK); // change to energy unit Ry when printing, S_tot, w/ model // wenfei add 2021/11/2 @@ -760,10 +764,10 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, if (!PARAM.inp.deepks_equiv) // training with stress label not supported by equivariant version now { - GlobalC::ld.cal_gvepsl(GlobalC::ucell.nat); + GlobalC::ld.cal_gvepsl(ucell.nat); LCAO_deepks_io::save_npy_gvepsl( - GlobalC::ucell.nat, + ucell.nat, GlobalC::ld.des_per_atom, GlobalC::ld.gvepsl_tensor, PARAM.globalv.global_out_dir, @@ -846,7 +850,8 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, // local pseudopotential, ewald, core correction, scc terms in force template -void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, +void Force_Stress_LCAO::calForcePwPart(const UnitCell& ucell, + ModuleBase::matrix& fvl_dvl, ModuleBase::matrix& fewalds, ModuleBase::matrix& fcc, ModuleBase::matrix& fscc, @@ -876,7 +881,7 @@ void Force_Stress_LCAO::calForcePwPart(ModuleBase::matrix& fvl_dvl, //-------------------------------------------------------- // force due to self-consistent charge. //-------------------------------------------------------- - f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, nlpp.numeric, GlobalC::ucell); + f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist, nlpp.numeric, ucell); return; } @@ -885,6 +890,7 @@ template <> void Force_Stress_LCAO::integral_part(const bool isGammaOnly, const bool isforce, const bool isstress, + const UnitCell& ucell, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -910,7 +916,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, flk.ftable(isforce, isstress, fsr, // mohan add 2024-06-15 - GlobalC::ucell, + ucell, psi, pelec, foverlap, @@ -935,6 +941,7 @@ template <> void Force_Stress_LCAO>::integral_part(const bool isGammaOnly, const bool isforce, const bool isstress, + const UnitCell& ucell, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi>* psi, @@ -959,7 +966,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn flk.ftable(isforce, isstress, fsr, // mohan add 2024-06-16 - GlobalC::ucell, + ucell, psi, pelec, foverlap, @@ -984,7 +991,8 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn // vlocal, hartree, ewald, core correction, exchange-correlation terms in stress template -void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, +void Force_Stress_LCAO::calStressPwPart(const UnitCell& ucell, + ModuleBase::matrix& sigmadvl, ModuleBase::matrix& sigmahar, ModuleBase::matrix& sigmaewa, ModuleBase::matrix& sigmacc, @@ -1022,7 +1030,7 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, //-------------------------------------------------------- for (int i = 0; i < 3; i++) { - sigmaxc(i, i) = -etxc / GlobalC::ucell.omega; + sigmaxc(i, i) = -etxc / ucell.omega; } // Exchange-correlation for PBE sc_pw.stress_gga(sigmaxc, rhopw, chr); @@ -1033,23 +1041,25 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, #include "module_base/mathzone.h" // do symmetry for total force template -void Force_Stress_LCAO::forceSymmetry(ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm) +void Force_Stress_LCAO::forceSymmetry(const UnitCell& ucell, + ModuleBase::matrix& fcs, + ModuleSymmetry::Symmetry* symm) { double d1, d2, d3; - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int iat = 0; iat < ucell.nat; iat++) { ModuleBase::Mathzone::Cartesian_to_Direct(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), - GlobalC::ucell.a1.x, - GlobalC::ucell.a1.y, - GlobalC::ucell.a1.z, - GlobalC::ucell.a2.x, - GlobalC::ucell.a2.y, - GlobalC::ucell.a2.z, - GlobalC::ucell.a3.x, - GlobalC::ucell.a3.y, - GlobalC::ucell.a3.z, + ucell.a1.x, + ucell.a1.y, + ucell.a1.z, + ucell.a2.x, + ucell.a2.y, + ucell.a2.z, + ucell.a3.x, + ucell.a3.y, + ucell.a3.z, d1, d2, d3); @@ -1059,20 +1069,20 @@ void Force_Stress_LCAO::forceSymmetry(ModuleBase::matrix& fcs, ModuleSymmetry fcs(iat, 2) = d3; } symm->symmetrize_vec3_nat(fcs.c); - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int iat = 0; iat < ucell.nat; iat++) { ModuleBase::Mathzone::Direct_to_Cartesian(fcs(iat, 0), fcs(iat, 1), fcs(iat, 2), - GlobalC::ucell.a1.x, - GlobalC::ucell.a1.y, - GlobalC::ucell.a1.z, - GlobalC::ucell.a2.x, - GlobalC::ucell.a2.y, - GlobalC::ucell.a2.z, - GlobalC::ucell.a3.x, - GlobalC::ucell.a3.y, - GlobalC::ucell.a3.z, + ucell.a1.x, + ucell.a1.y, + ucell.a1.z, + ucell.a2.x, + ucell.a2.y, + ucell.a2.z, + ucell.a3.x, + ucell.a3.y, + ucell.a3.z, d1, d2, d3); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index 35a596bdcb..927d08790e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -33,6 +33,7 @@ class Force_Stress_LCAO const bool isstress, const bool istestf, const bool istests, + const UnitCell& ucell, Parallel_Orbitals& pv, const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -59,9 +60,12 @@ class Force_Stress_LCAO Stress_Func sc_pw; Forces f_pw; - void forceSymmetry(ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm); + void forceSymmetry(const UnitCell& ucell, + ModuleBase::matrix& fcs, + ModuleSymmetry::Symmetry* symm); - void calForcePwPart(ModuleBase::matrix& fvl_dvl, + void calForcePwPart(const UnitCell& ucell, + ModuleBase::matrix& fvl_dvl, ModuleBase::matrix& fewalds, ModuleBase::matrix& fcc, ModuleBase::matrix& fscc, @@ -76,6 +80,7 @@ class Force_Stress_LCAO void integral_part(const bool isGammaOnly, const bool isforce, const bool isstress, + const UnitCell& ucell, ForceStressArrays& fsr, // mohan add 2024-06-15 const elecstate::ElecState* pelec, const psi::Psi* psi, @@ -97,7 +102,8 @@ class Force_Stress_LCAO const Parallel_Orbitals& pv, const K_Vectors& kv); - void calStressPwPart(ModuleBase::matrix& sigmadvl, + void calStressPwPart(const UnitCell& ucell, + ModuleBase::matrix& sigmadvl, ModuleBase::matrix& sigmahar, ModuleBase::matrix& sigmaewa, ModuleBase::matrix& sigmacc, From 3116e748ae6e4344d7025eb2e5eb001ffeca8be1 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 11:39:41 +0800 Subject: [PATCH 03/10] change GlobalC::ucell in the grid_init.cpp --- source/module_esolver/lcao_before_scf.cpp | 2 +- source/module_esolver/lcao_others.cpp | 2 +- source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h | 1 + source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp | 6 +++--- 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index b7ec3ce705..e7b7e5f6ac 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -174,7 +174,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) } // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big); + LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big); // init Hamiltonian if (this->p_hamilt != nullptr) diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 10d688dc00..6518717a6f 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -175,7 +175,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, orb_, *this->pw_rho, *this->pw_big); + LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell,orb_, *this->pw_rho, *this->pw_big); // init Hamiltonian if (this->p_hamilt != nullptr) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h index f98f2433ab..f1e4979642 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h @@ -41,6 +41,7 @@ void build_Nonlocal_mu_new(const Parallel_Orbitals& pv, void grid_prepare(const Grid_Technique& gt, Gint_Gamma& gint_gamma, Gint_k& gint_k, + const UnitCell& ucell, const LCAO_Orbitals& orb, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp index 0a1997c9d2..d14e3ef8a9 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/grid_init.cpp @@ -15,13 +15,13 @@ void grid_prepare( const Grid_Technique& gt, Gint_Gamma &gint_gamma, Gint_k &gint_k, + const UnitCell& ucell, const LCAO_Orbitals& orb, const ModulePW::PW_Basis& rhopw, const ModulePW::PW_Basis_Big& bigpw) { ModuleBase::TITLE("LCAO_domain","grid_prepare"); ModuleBase::timer::tick("LCAO_domain","grid_prepare"); - const UnitCell* ucell = &GlobalC::ucell; if(PARAM.globalv.gamma_only_local) { gint_gamma.prep_grid( @@ -39,7 +39,7 @@ void grid_prepare( rhopw.ny, rhopw.nplane, rhopw.startz_current, - ucell, + &ucell, &orb); } else // multiple k-points @@ -60,7 +60,7 @@ void grid_prepare( rhopw.ny, rhopw.nplane, rhopw.startz_current, - ucell, + &ucell, &orb); } From 79f15d588868292e6aaca50da58b1987f09775a9 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 11:46:38 +0800 Subject: [PATCH 04/10] update ucell in hamilt_lcao --- source/module_esolver/esolver_gets.cpp | 6 ++- source/module_esolver/lcao_before_scf.cpp | 1 + source/module_esolver/lcao_others.cpp | 1 + .../hamilt_lcaodft/hamilt_lcao.cpp | 39 +++++++++++-------- .../hamilt_lcaodft/hamilt_lcao.h | 7 +++- 5 files changed, 34 insertions(+), 20 deletions(-) diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 7cd0258621..3aed24d99d 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -108,7 +108,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) if (PARAM.inp.nspin == 4) { this->p_hamilt - = new hamilt::HamiltLCAO, std::complex>(&this->pv, + = new hamilt::HamiltLCAO, std::complex>(ucell, + &this->pv, this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); @@ -117,7 +118,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) } else { - this->p_hamilt = new hamilt::HamiltLCAO, double>(&this->pv, + this->p_hamilt = new hamilt::HamiltLCAO, double>(ucell, + &this->pv, this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index e7b7e5f6ac..213de30a60 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -188,6 +188,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->p_hamilt = new hamilt::HamiltLCAO( PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + ucell, &this->pv, this->pelec->pot, this->kv, diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 6518717a6f..6194e20284 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -189,6 +189,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->p_hamilt = new hamilt::HamiltLCAO( PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + ucell, &this->pv, this->pelec->pot, this->kv, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index dfc9d0a55a..ca76150f17 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -43,7 +43,11 @@ namespace hamilt { template -HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff) +HamiltLCAO::HamiltLCAO(const UnitCell& ucell, + const Parallel_Orbitals* paraV, + const K_Vectors& kv_in, + const TwoCenterIntegrator& intor_overlap_orb, + const std::vector& orb_cutoff) { this->classname = "HamiltLCAO"; @@ -58,7 +62,7 @@ HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& this->kv->kvec_d, this->hR, this->sR, - &GlobalC::ucell, + &ucell, orb_cutoff, &GlobalC::GridD, &intor_overlap_orb); @@ -67,6 +71,7 @@ HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& template HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Gint_k* GK_in, + const UnitCell& ucell, const Parallel_Orbitals* paraV, elecstate::Potential* pot_in, const K_Vectors& kv_in, @@ -133,7 +138,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, this->sR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); @@ -144,7 +149,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* ekinetic = new EkineticNew>(this->hsk, this->kv->kvec_d, this->hR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.kinetic_orb.get()); @@ -158,7 +163,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* nonlocal = new NonlocalNew>(this->hsk, this->kv->kvec_d, this->hR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get()); @@ -180,7 +185,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, pot_in, this->hR, // no explicit call yet - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, PARAM.inp.nspin @@ -195,7 +200,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* deepks = new DeePKS>(this->hsk, this->kv->kvec_d, this->hR, // no explicit call yet - &GlobalC::ucell, + &ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_alpha.get(), &orb, @@ -221,7 +226,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, dftu = new DFTU>(this->hsk, this->kv->kvec_d, this->hR, - GlobalC::ucell, + ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), @@ -249,7 +254,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, kv->kvec_d, pot_in, this->hR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, PARAM.inp.nspin); @@ -265,7 +270,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, this->sR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); @@ -285,7 +290,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* ekinetic = new EkineticNew>(this->hsk, this->kv->kvec_d, this->hR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.kinetic_orb.get()); @@ -299,7 +304,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* nonlocal = new NonlocalNew>(this->hsk, this->kv->kvec_d, this->hR, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb_beta.get()); @@ -321,7 +326,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* deepks = new DeePKS>(this->hsk, this->kv->kvec_d, hR, - &GlobalC::ucell, + &ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_alpha.get(), &orb, @@ -338,7 +343,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* td_ekinetic = new TDEkinetic>(this->hsk, this->hR, kv, - &GlobalC::ucell, + &ucell, orb.cutoffs(), &GlobalC::GridD, two_center_bundle.overlap_orb.get()); @@ -347,7 +352,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* td_nonlocal = new TDNonlocal>(this->hsk, this->kv->kvec_d, this->hR, - &GlobalC::ucell, + &ucell, orb, &GlobalC::GridD); this->getOperator()->add(td_nonlocal); @@ -367,7 +372,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, dftu = new DFTU>(this->hsk, this->kv->kvec_d, this->hR, - GlobalC::ucell, + ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), @@ -382,7 +387,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hsk, kv->kvec_d, this->hR, - GlobalC::ucell, + ucell, &GlobalC::GridD, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs() diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 0cd597b5b0..0c8a5bc188 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -30,6 +30,7 @@ class HamiltLCAO : public Hamilt using TAC = std::pair>; HamiltLCAO(Gint_Gamma* GG_in, Gint_k* GK_in, + const UnitCell& ucell, const Parallel_Orbitals* paraV, elecstate::Potential* pot_in, const K_Vectors& kv_in, @@ -46,7 +47,11 @@ class HamiltLCAO : public Hamilt /** * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer */ - HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff); + HamiltLCAO(const UnitCell& ucell, + const Parallel_Orbitals* paraV, + const K_Vectors& kv_in, + const TwoCenterIntegrator& intor_overlap_orb, + const std::vector& orb_cutoff); ~HamiltLCAO() { From b0e03525b5dc60d515646da5ab42eb86df1a7f61 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 15:30:27 +0800 Subject: [PATCH 05/10] update ucell in LCAO_alloacte --- source/module_esolver/esolver_ks_lcao.cpp | 2 +- .../hamilt_lcaodft/LCAO_allocate.cpp | 16 ++++++++-------- .../hamilt_lcaodft/LCAO_domain.h | 2 +- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 9331c206ff..b5a62ac840 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -164,7 +164,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // 5) initialize Hamilt in LCAO // * allocate H and S matrices according to computational resources // * set the 'trace' between local H/S and global H/S - LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, pv, this->kv.get_nks(), orb_); + LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell , pv, this->kv.get_nks(), orb_); #ifdef __EXX // 6) initialize exx diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp index b9e0be5a46..2b9969554e 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp @@ -7,7 +7,7 @@ namespace LCAO_domain { -void divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv,const int& nks, const LCAO_Orbitals& orb) { +void divide_HS_in_frag(const bool isGamma, const UnitCell& ucell, Parallel_Orbitals& pv,const int& nks, const LCAO_Orbitals& orb) { ModuleBase::TITLE("LCAO_domain", "divide_HS_in_frag"); //(1), (2): set up matrix division have been moved into ESolver_KS_LCAO::init_basis_lcao @@ -19,22 +19,22 @@ void divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv,const int& nks, if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) { // allocate relevant data structures for calculating descriptors std::vector na; - na.resize(GlobalC::ucell.ntype); - for (int it = 0; it < GlobalC::ucell.ntype; it++) { - na[it] = GlobalC::ucell.atoms[it].na; + na.resize(ucell.ntype); + for (int it = 0; it < ucell.ntype; it++) { + na[it] = ucell.atoms[it].na; } GlobalC::ld.init(orb, - GlobalC::ucell.nat, - GlobalC::ucell.ntype, + ucell.nat, + ucell.ntype, pv, na); if (PARAM.inp.deepks_scf) { if (isGamma) { - GlobalC::ld.allocate_V_delta(GlobalC::ucell.nat); + GlobalC::ld.allocate_V_delta(ucell.nat); } else { - GlobalC::ld.allocate_V_delta(GlobalC::ucell.nat, nks); + GlobalC::ld.allocate_V_delta(ucell.nat, nks); } } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h index f1e4979642..f7a581b16c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h @@ -174,7 +174,7 @@ void build_ST_new(ForceStressArrays& fsr, */ void zeros_HSR(const char& mtype, LCAO_HS_Arrays& HS_arrays); -void divide_HS_in_frag(const bool isGamma, Parallel_Orbitals& pv, const int& nks, const LCAO_Orbitals& orb); +void divide_HS_in_frag(const bool isGamma, const UnitCell& ucell, Parallel_Orbitals& pv, const int& nks, const LCAO_Orbitals& orb); template void set_mat2d(const int& global_ir, From 5b77d7f0e926cced2a579551abaa0659bfd72a0d Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 15:46:13 +0800 Subject: [PATCH 06/10] update ucell in record_adj --- source/module_esolver/esolver_gets.cpp | 2 +- source/module_esolver/lcao_before_scf.cpp | 2 +- source/module_esolver/lcao_others.cpp | 2 +- .../hamilt_lcaodft/record_adj.cpp | 126 +++++++++--------- .../hamilt_lcaodft/record_adj.h | 4 +- 5 files changed, 68 insertions(+), 68 deletions(-) diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 3aed24d99d..3571d4788f 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -101,7 +101,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) PARAM.inp.test_atom_input); Record_adj RA; - RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); if (this->p_hamilt == nullptr) { diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 213de30a60..3d018753de 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -127,7 +127,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 6194e20284..43753ce01b 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -128,7 +128,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - this->RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp index 43569e2e24..695e6327db 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp @@ -38,20 +38,20 @@ void Record_adj::delete_grid() // If multi-k, calculate nnr at the same time. // be called only once in an ion-step. //-------------------------------------------- -void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff) +void Record_adj::for_2d(const UnitCell& ucell, Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff) { ModuleBase::TITLE("Record_adj","for_2d"); ModuleBase::timer::tick("Record_adj","for_2d"); - assert(GlobalC::ucell.nat > 0); + assert(ucell.nat > 0); if (!gamma_only) { delete[] pv.nlocdim; delete[] pv.nlocstart; - pv.nlocdim = new int[GlobalC::ucell.nat]; - pv.nlocstart = new int[GlobalC::ucell.nat]; - ModuleBase::GlobalFunc::ZEROS(pv.nlocdim, GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pv.nlocstart, GlobalC::ucell.nat); + pv.nlocdim = new int[ucell.nat]; + pv.nlocstart = new int[ucell.nat]; + ModuleBase::GlobalFunc::ZEROS(pv.nlocdim, ucell.nat); + ModuleBase::GlobalFunc::ZEROS(pv.nlocstart, ucell.nat); pv.nnr = 0; } { @@ -59,7 +59,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto ModuleBase::Vector3 tau1, tau2, dtau; ModuleBase::Vector3 dtau1, dtau2, tau0; - this->na_proc = GlobalC::ucell.nat; + this->na_proc = ucell.nat; // number of adjacents for each atom. this->na_each = new int[na_proc]; @@ -67,15 +67,15 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto int iat = 0; - for (int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1) + for (int T1 = 0; T1 < ucell.ntype; ++T1) { - Atom* atom1 = &GlobalC::ucell.atoms[T1]; + Atom* atom1 = &ucell.atoms[T1]; for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; //GlobalC::GridD.Find_atom( tau1 ); - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1 ,T1, I1); - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + GlobalC::GridD.Find_atom(ucell, tau1 ,T1, I1); + const int start1 = ucell.itiaiw2iwt(T1, I1, 0); if(!gamma_only) { pv.nlocstart[iat] = pv.nnr; } @@ -84,10 +84,10 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto { const int T2 = GlobalC::GridD.getType(ad); const int I2 = GlobalC::GridD.getNatom(ad); - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + const int start2 = ucell.itiaiw2iwt(T2, I2, 0); tau2 = GlobalC::GridD.getAdjacentTau(ad); dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff[T1] + orb_cutoff[T2]; bool is_adj = false; @@ -102,17 +102,17 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto { const int T0 = GlobalC::GridD.getType(ad0); //const int I0 = GlobalC::GridD.getNatom(ad0); - //const int iat0 = GlobalC::ucell.itia2iat(T0, I0); - //const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + //const int iat0 = ucell.itia2iat(T0, I0); + //const int start0 = ucell.itiaiw2iwt(T0, I0, 0); tau0 = GlobalC::GridD.getAdjacentTau(ad0); dtau1 = tau0 - tau1; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double distance1 = dtau1.norm() * ucell.lat0; + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); dtau2 = tau0 - tau2; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double distance2 = dtau2.norm() * ucell.lat0; + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); if( distance1 < rcut1 && distance2 < rcut2 ) { @@ -135,7 +135,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto if(mu<0) {continue; } - for(int jj=0; jjtau[I1]; //GlobalC::GridD.Find_atom( tau1 ); AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1 ,T1, I1, &adjs); + GlobalC::GridD.Find_atom(ucell, tau1 ,T1, I1, &adjs); // (2) search among all adjacent atoms. int cb = 0; @@ -209,7 +209,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto const int I2 = adjs.natom[ad]; tau2 = adjs.adjacent_tau[ad]; dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff[T1] + orb_cutoff[T2]; @@ -221,17 +221,17 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto { const int T0 = adjs.ntype[ad0]; //const int I0 = GlobalC::GridD.getNatom(ad0); - //const int iat0 = GlobalC::ucell.itia2iat(T0, I0); - //const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + //const int iat0 = ucell.itia2iat(T0, I0); + //const int start0 = ucell.itiaiw2iwt(T0, I0, 0); tau0 = adjs.adjacent_tau[ad0]; dtau1 = tau0 - tau1; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double distance1 = dtau1.norm() * ucell.lat0; + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); dtau2 = tau0 - tau2; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double distance2 = dtau2.norm() * ucell.lat0; + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); if( distance1 < rcut1 && distance2 < rcut2 ) { @@ -267,14 +267,14 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto // This will record the orbitals according to // grid division (cut along z direction) //-------------------------------------------- -void Record_adj::for_grid(const Grid_Technique >, const std::vector& orb_cutoff) +void Record_adj::for_grid(const UnitCell& ucell, const Grid_Technique >, const std::vector& orb_cutoff) { ModuleBase::TITLE("Record_adj","for_grid"); ModuleBase::timer::tick("Record_adj","for_grid"); this->na_proc = 0; - this->iat2ca = new int[GlobalC::ucell.nat]; - for(int iat=0; iatiat2ca = new int[ucell.nat]; + for(int iat=0; iat& o #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int iat=0; iat& o tau1 = atom1->tau[I1]; //GlobalC::GridD.Find_atom(tau1); AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); + GlobalC::GridD.Find_atom(ucell, tau1, T1, I1, &adjs); for (int ad = 0; ad < adjs.adj_num+1; ad++) { const int T2 = adjs.ntype[ad]; const int I2 = adjs.natom[ad]; - const int iat2 = GlobalC::ucell.itia2iat(T2, I2); + const int iat2 = ucell.itia2iat(T2, I2); if(gt.in_this_processor[iat2]) { - //Atom* atom2 = &GlobalC::ucell.atoms[T2]; + //Atom* atom2 = &ucell.atoms[T2]; tau2 = adjs.adjacent_tau[ad]; dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff[T1] + orb_cutoff[T2]; @@ -339,18 +339,18 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o { const int T0 = GlobalC::GridD.getType(ad0); const int I0 = GlobalC::GridD.getNatom(ad0); - const int iat0 = GlobalC::ucell.itia2iat(T0, I0); - const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + const int iat0 = ucell.itia2iat(T0, I0); + const int start0 = ucell.itiaiw2iwt(T0, I0, 0); tau0 = GlobalC::GridD.getAdjacentTau(ad0); dtau1 = tau0 - tau1; dtau2 = tau0 - tau2; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; + double distance1 = dtau1.norm() * ucell.lat0; + double distance2 = dtau2.norm() * ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); if( distance1 < rcut1 && distance2 < rcut2 ) { @@ -390,11 +390,11 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int iat=0; iat& o tau1 = atom1->tau[I1]; //GlobalC::GridD.Find_atom(tau1); AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); + GlobalC::GridD.Find_atom(ucell, tau1, T1, I1, &adjs); int cb = 0; for (int ad = 0; ad < adjs.adj_num+1; ad++) { const int T2 = adjs.ntype[ad]; const int I2 = adjs.natom[ad]; - const int iat2 = GlobalC::ucell.itia2iat(T2, I2); + const int iat2 = ucell.itia2iat(T2, I2); // key of this function if(gt.in_this_processor[iat2]) { - //Atom* atom2 = &GlobalC::ucell.atoms[T2]; + //Atom* atom2 = &ucell.atoms[T2]; tau2 = adjs.adjacent_tau[ad]; dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff[T1] + orb_cutoff[T2]; // check the distance @@ -439,18 +439,18 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o { const int T0 = GlobalC::GridD.getType(ad0); const int I0 = GlobalC::GridD.getNatom(ad0); - const int iat0 = GlobalC::ucell.itia2iat(T0, I0); - const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + const int iat0 = ucell.itia2iat(T0, I0); + const int start0 = ucell.itiaiw2iwt(T0, I0, 0); tau0 = GlobalC::GridD.getAdjacentTau(ad0); dtau1 = tau0 - tau1; dtau2 = tau0 - tau2; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; + double distance1 = dtau1.norm() * ucell.lat0; + double distance2 = dtau2.norm() * ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); if( distance1 < rcut1 && distance2 < rcut2 ) { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h index d2f214632f..06cf238a81 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h @@ -21,13 +21,13 @@ class Record_adj // This will record the orbitals according to // HPSEPS's 2D block division. //-------------------------------------------- - void for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff); + void for_2d(const UnitCell& ucell, Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff); //-------------------------------------------- // This will record the orbitals according to // grid division (cut along z direction) //-------------------------------------------- - void for_grid(const Grid_Technique >, const std::vector& orb_cutoff); + void for_grid(const UnitCell& ucell, const Grid_Technique >, const std::vector& orb_cutoff); void delete_grid(); From b838e994b0532805620eadafd994667e44c5b718 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 15:56:35 +0800 Subject: [PATCH 07/10] update ucell in spra_dh.cpp --- .../hamilt_lcaodft/spar_dh.cpp | 38 ++++++++++--------- .../hamilt_lcaodft/spar_dh.h | 6 ++- .../module_hamilt_lcao/module_gint/gint_k.h | 2 +- .../module_gint/gint_k_sparse1.cpp | 2 +- source/module_io/output_mat_sparse.cpp | 1 + source/module_io/write_HS_R.cpp | 7 +++- source/module_io/write_HS_R.h | 1 + 7 files changed, 33 insertions(+), 24 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp index 3487276c5f..209f067241 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.cpp @@ -4,7 +4,8 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include -void sparse_format::cal_dH(const Parallel_Orbitals& pv, +void sparse_format::cal_dH(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const TwoCenterBundle& two_center_bundle, @@ -36,7 +37,7 @@ void sparse_format::cal_dH(const Parallel_Orbitals& pv, 'T', cal_deri, cal_stress, - GlobalC::ucell, + ucell, orb, pv, two_center_bundle, @@ -48,19 +49,19 @@ void sparse_format::cal_dH(const Parallel_Orbitals& pv, fsr_dh, nullptr, true, - GlobalC::ucell, + ucell, orb, *(two_center_bundle.overlap_orb_beta), &GlobalC::GridD); - sparse_format::cal_dSTN_R(pv, HS_Arrays, fsr_dh, grid, orb.cutoffs(), current_spin, sparse_thr); + sparse_format::cal_dSTN_R(ucell,pv, HS_Arrays, fsr_dh, grid, orb.cutoffs(), current_spin, sparse_thr); delete[] fsr_dh.DHloc_fixedR_x; delete[] fsr_dh.DHloc_fixedR_y; delete[] fsr_dh.DHloc_fixedR_z; gint_k - .cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, GlobalC::ucell, GlobalC::GridD); + .cal_dvlocal_R_sparseMatrix(current_spin, sparse_thr, HS_Arrays, &pv, ucell, GlobalC::GridD); return; } @@ -90,7 +91,8 @@ void sparse_format::set_R_range(std::set>& all_R_coor, return; } -void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, +void sparse_format::cal_dSTN_R(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, ForceStressArrays& fsr, Grid_Driver& grid, @@ -107,25 +109,25 @@ void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, double temp_value_double; std::complex temp_value_complex; - for (int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1) + for (int T1 = 0; T1 < ucell.ntype; ++T1) { - Atom* atom1 = &GlobalC::ucell.atoms[T1]; + Atom* atom1 = &ucell.atoms[T1]; for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; - grid.Find_atom(GlobalC::ucell, tau1, T1, I1); - Atom* atom1 = &GlobalC::ucell.atoms[T1]; - const int start = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + grid.Find_atom(ucell, tau1, T1, I1); + Atom* atom1 = &ucell.atoms[T1]; + const int start = ucell.itiaiw2iwt(T1, I1, 0); for (int ad = 0; ad < grid.getAdjacentNum() + 1; ++ad) { const int T2 = grid.getType(ad); const int I2 = grid.getNatom(ad); - Atom* atom2 = &GlobalC::ucell.atoms[T2]; + Atom* atom2 = &ucell.atoms[T2]; tau2 = grid.getAdjacentTau(ad); dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = orb_cutoff[T1] + orb_cutoff[T2]; bool adj = false; @@ -144,11 +146,11 @@ void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, dtau1 = tau0 - tau1; dtau2 = tau0 - tau2; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; + double distance1 = dtau1.norm() * ucell.lat0; + double distance2 = dtau2.norm() * ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); if (distance1 < rcut1 && distance2 < rcut2) { @@ -160,7 +162,7 @@ void sparse_format::cal_dSTN_R(const Parallel_Orbitals& pv, if (adj) { - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + const int start2 = ucell.itiaiw2iwt(T2, I2, 0); Abfs::Vector3_Order dR(grid.getBox(ad).x, grid.getBox(ad).y, grid.getBox(ad).z); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h index ce855a53dc..f0657dfdd7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_dh.h @@ -11,7 +11,8 @@ namespace sparse_format { -void cal_dH(const Parallel_Orbitals& pv, +void cal_dH(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, const TwoCenterBundle& two_center_bundle, @@ -24,7 +25,8 @@ void cal_dH(const Parallel_Orbitals& pv, void set_R_range(std::set>& all_R_coor, Grid_Driver& grid); // be called by 'cal_dH_sparse' -void cal_dSTN_R(const Parallel_Orbitals& pv, +void cal_dSTN_R(const UnitCell& ucell, + const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, ForceStressArrays& fsr, // mohan add 2024-06-16 Grid_Driver& grid, diff --git a/source/module_hamilt_lcao/module_gint/gint_k.h b/source/module_hamilt_lcao/module_gint/gint_k.h index 3aa368e286..21968b3a27 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k.h +++ b/source/module_hamilt_lcao/module_gint/gint_k.h @@ -95,7 +95,7 @@ class Gint_k : public Gint { const double& sparse_threshold, LCAO_HS_Arrays& HS_Arrays, const Parallel_Orbitals* pv, - UnitCell& ucell, + const UnitCell& ucell, Grid_Driver& gdriver); private: diff --git a/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp b/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp index 952ed5c02b..753c1eff17 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_sparse1.cpp @@ -322,7 +322,7 @@ void Gint_k::cal_dvlocal_R_sparseMatrix(const int& current_spin, const double& sparse_threshold, LCAO_HS_Arrays& HS_Arrays, const Parallel_Orbitals* pv, - UnitCell& ucell, + const UnitCell& ucell, Grid_Driver& gdriver) { ModuleBase::TITLE("Gint_k", "cal_dvlocal_R_sparseMatrix"); diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index 9891e06836..b5139876a5 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -60,6 +60,7 @@ void output_mat_sparse(const bool& out_mat_hsR, output_dHR(istep, v_eff, gint_k, // mohan add 2024-04-01 + ucell, pv, HS_Arrays, grid, // mohan add 2024-04-06 diff --git a/source/module_io/write_HS_R.cpp b/source/module_io/write_HS_R.cpp index 590f95fb78..5b59f2ee24 100644 --- a/source/module_io/write_HS_R.cpp +++ b/source/module_io/write_HS_R.cpp @@ -86,6 +86,7 @@ void ModuleIO::output_HSR(const int& istep, void ModuleIO::output_dHR(const int& istep, const ModuleBase::matrix& v_eff, Gint_k& gint_k, // mohan add 2024-04-01 + const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, // mohan add 2024-04-06 @@ -105,7 +106,8 @@ void ModuleIO::output_dHR(const int& istep, // mohan add 2024-04-01 const int cspin = 0; - sparse_format::cal_dH(pv, + sparse_format::cal_dH(ucell, + pv, HS_Arrays, grid, two_center_bundle, @@ -129,7 +131,8 @@ void ModuleIO::output_dHR(const int& istep, } } - sparse_format::cal_dH(pv, + sparse_format::cal_dH(ucell, + pv, HS_Arrays, grid, two_center_bundle, diff --git a/source/module_io/write_HS_R.h b/source/module_io/write_HS_R.h index 5a23a6796e..6c9bd46768 100644 --- a/source/module_io/write_HS_R.h +++ b/source/module_io/write_HS_R.h @@ -31,6 +31,7 @@ void output_HSR(const int& istep, void output_dHR(const int& istep, const ModuleBase::matrix& v_eff, Gint_k& gint_k, // mohan add 2024-04-01 + const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, Grid_Driver& grid, // mohan add 2024-04-06 From d5e8b1398e7906efd3adbf8e52d3b55f6a6edb94 Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 16:40:34 +0800 Subject: [PATCH 08/10] update ucell in wavefunc_in_pw --- .../hamilt_lcaodft/wavefunc_in_pw.cpp | 79 ++++++++++--------- .../hamilt_lcaodft/wavefunc_in_pw.h | 5 +- 2 files changed, 47 insertions(+), 37 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp index db06c13f3a..9faf6437e8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp @@ -9,17 +9,18 @@ #include "module_hamilt_pw/hamilt_pwdft/soc.h" void Wavefunc_in_pw::make_table_q( + const UnitCell &ucell, std::vector &fn, ModuleBase::realArray &table_local) { ModuleBase::TITLE("Wavefunc_in_pw","make_table_q"); - if( fn.size() != static_cast(GlobalC::ucell.ntype) ) + if( fn.size() != static_cast(ucell.ntype) ) { ModuleBase::WARNING_QUIT("Wavefunc_in_pw::make_table_q","maybe NUMERICAL_ORBITAL is not read in, please check."); } - for(int it=0; it=0); const int npw = wfc_basis->npwk[ik]; - const int total_lm = ( GlobalC::ucell.lmax + 1) * ( GlobalC::ucell.lmax + 1); + const int total_lm = ( ucell.lmax + 1) * ( ucell.lmax + 1); ModuleBase::matrix ylm(total_lm, npw); std::complex *aux = new std::complex[npw]; double *chiaux = nullptr; @@ -273,23 +280,23 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, //int index = 0; double *flq = new double[npw]; int iwall=0; - for (int it = 0;it < GlobalC::ucell.ntype;it++) + for (int it = 0;it < ucell.ntype;it++) { - for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++) + for (int ia = 0;ia < ucell.atoms[it].na;ia++) { std::complex* sk = sf.get_sk(ik, it, ia, wfc_basis); int ic = 0; - for(int L = 0; L < GlobalC::ucell.atoms[it].nwl+1; L++) + for(int L = 0; L < ucell.atoms[it].nwl+1; L++) { std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, L); //mohan 2010-04-19 - for(int N=0; N < GlobalC::ucell.atoms[it].l_nchi[L]; N++) + for(int N=0; N < ucell.atoms[it].l_nchi[L]; N++) { // GlobalV::ofs_running << " it=" << it << " ia=" << ia << " L=" << L << " N=" << N << std::endl; for(int ig=0; ig GlobalC::ucell.natomwfc) + if (iwall + 2 * L + 1 > ucell.natomwfc) ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs"); for (int ig = 0; ig < npw; ig++) { @@ -389,19 +396,19 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, } iwall += 2*L +1; } // end else INPUT.starting_spin_angle || !PARAM.globalv.domag - } // end if GlobalC::ucell.atoms[it].has_so + } // end if ucell.atoms[it].has_so else {//atomic_wfc_nc double alpha, gamman; std::complex fup, fdown; - //alpha = GlobalC::ucell.magnet.angle1_[it]; - //gamman = -GlobalC::ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; - alpha = GlobalC::ucell.atoms[it].angle1[ia]; - gamman = -GlobalC::ucell.atoms[it].angle2[ia] + 0.5*ModuleBase::PI; + //alpha = ucell.magnet.angle1_[it]; + //gamman = -ucell.magnet.angle2_[it] + 0.5*ModuleBase::PI; + alpha = ucell.atoms[it].angle1[ia]; + gamman = -ucell.atoms[it].angle2[ia] + 0.5*ModuleBase::PI; for(int m = 0;m<2*L+1;m++) { const int lm = L*L +m; - if (iwall + 2 * L + 1 > GlobalC::ucell.natomwfc) + if (iwall + 2 * L + 1 > ucell.natomwfc) { ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs"); } @@ -430,7 +437,7 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const int& ik, iwall++; } // end m iwall += 2*L+1; - } // end else GlobalC::ucell.atoms[it].has_so + } // end else ucell.atoms[it].has_so } // end for is_N } // end if PARAM.inp.noncolin else diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h index e716483b57..47cd5f9cff 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.h @@ -17,10 +17,12 @@ namespace Wavefunc_in_pw { void make_table_q( + const UnitCell &ucell, std::vector &orbital_files, ModuleBase::realArray &table_local); void integral( + const UnitCell& ucell, const int meshr, // number of mesh points const double *psir, const double *r, @@ -34,7 +36,8 @@ namespace Wavefunc_in_pw const double &ecut, const double &beta); - void produce_local_basis_in_pw(const int& ik, + void produce_local_basis_in_pw(const UnitCell& ucell, + const int& ik, const ModulePW::PW_Basis_K* wfc_basis, const Structure_Factor& sf, ModuleBase::ComplexMatrix& psi, From 519c09fa23e69d7f4f8a2111343cf33d46ce405d Mon Sep 17 00:00:00 2001 From: A-006 <3158793232@qq.com> Date: Fri, 29 Nov 2024 16:47:42 +0800 Subject: [PATCH 09/10] fix bug in wavefunc --- source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp index 9faf6437e8..9e9d49d433 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp @@ -134,7 +134,7 @@ void Wavefunc_in_pw::make_table_q( } } double* table = new double[PARAM.globalv.nqx]; - Wavefunc_in_pw::integral(meshr, psir, radial, rab, L, table); + Wavefunc_in_pw::integral(ucell,meshr, psir, radial, rab, L, table); for(int iq=0; iq Date: Fri, 29 Nov 2024 10:27:29 +0000 Subject: [PATCH 10/10] [pre-commit.ci lite] apply automatic fixes --- .../hamilt_lcaodft/wavefunc_in_pw.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp index 9e9d49d433..b26e85cf1f 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp @@ -304,7 +304,8 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const UnitCell& ucell, /* for(int is_N = 0; is_N < 2; is_N++)*/ //for rotate base for(int is_N = 0; is_N < 1; is_N++) { - if(L==0 && is_N==1) continue; + if(L==0 && is_N==1) { continue; +} if(ucell.atoms[it].ncpp.has_so) { const double j = std::abs(double(L+is_N) - 0.5); @@ -332,15 +333,16 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const UnitCell& ucell, std::complex fup,fdown; //int nc; //This routine creates two functions only in the case j=l+1/2 or exit in the other case - if(fabs(j-L+0.5)<1e-4) continue; + if(fabs(j-L+0.5)<1e-4) { continue; +} delete[] chiaux; chiaux = new double [npw]; //Find the functions j= l- 1/2 - if(L==0) + if(L==0) { for(int ig=0;ig ucell.natomwfc) + if (iwall + 2 * L + 1 > ucell.natomwfc) { ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs"); +} for (int ig = 0; ig < npw; ig++) { aux[ig] = sk[ig] * ylm(lm,ig) * chiaux[ig];