diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 5a39566c49..8c947a9fb2 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -136,8 +136,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) if (PARAM.inp.out_mat_r) { cal_r_overlap_R r_matrix; - r_matrix.init(pv, orb_); - r_matrix.out_rR(istep); + r_matrix.init(ucell,pv, orb_); + r_matrix.out_rR(ucell,istep); } ModuleBase::timer::tick("ESolver_GetS", "runner"); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 31d08fb71f..67507d7831 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -991,6 +991,8 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.out_dm1, false, PARAM.inp.out_app_flag, + ucell.get_iat2iwt(), + &ucell.nat, istep); //! 7) write density matrix @@ -1232,7 +1234,8 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.nnkpfile, PARAM.inp.wannier_spin); - myWannier.calculate(this->pelec->ekb, + myWannier.calculate(ucell, + this->pelec->ekb, this->pw_wfc, this->pw_big, this->sf, @@ -1251,7 +1254,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) PARAM.inp.wannier_spin, orb_); - myWannier.calculate(this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); + myWannier.calculate(ucell,this->pelec->ekb, this->kv, *(this->psi), &(this->pv)); } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90"); } @@ -1263,12 +1266,13 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation"); berryphase bp(&(this->pv)); - bp.lcao_init(this->kv, + bp.lcao_init(ucell, + this->kv, this->GridT, orb_); // additional step before calling // macroscopic_polarization (why capitalize // the function name?) - bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation"); } } diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index ef231e84e9..5a6d4c4087 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -349,7 +349,7 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) { std::stringstream ss_dipole; ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; - ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); + ModuleIO::write_dipole(ucell,pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); } } if (TD_Velocity::out_current == true) @@ -357,7 +357,8 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) elecstate::DensityMatrix, double>* tmp_DM = dynamic_cast>*>(this->pelec)->get_DM(); - ModuleIO::write_current(istep, + ModuleIO::write_current(ucell, + istep, this->psi, pelec, kv, diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 67e4d33dd1..c0e91415f4 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -585,8 +585,8 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) PARAM.inp.out_wannier_wvfn_formatted, PARAM.inp.nnkpfile, PARAM.inp.wannier_spin); - - wan.calculate(this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi); + wan.set_tpiba_omega(ucell.tpiba, ucell.omega); + wan.calculate(ucell,this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation"); } @@ -595,7 +595,7 @@ void ESolver_KS_PW::after_scf(UnitCell& ucell, const int istep) { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase polarization"); berryphase bp; - bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); + bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv); std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase polarization"); } } @@ -783,7 +783,7 @@ void ESolver_KS_PW::after_all_runners(UnitCell& ucell) //! 6) Print out electronic wave functions in real space if (PARAM.inp.out_wfc_r == 1) // Peize Lin add 2021.11.21 { - ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv); + ModuleIO::write_psi_r_1(ucell,this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv); } //! 7) Use Kubo-Greenwood method to compute conductivities diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index c5c2aae021..50012599c1 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -48,7 +48,10 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) if (cal_type == "test_memory") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "testing memory"); - Cal_Test::test_memory(this->pw_rho, + Cal_Test::test_memory(ucell.nat, + ucell.ntype, + ucell.GGT, + this->pw_rho, this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); @@ -324,7 +327,8 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) IState_Envelope IEP(this->pelec); if (PARAM.globalv.gamma_only_local) { - IEP.begin(this->psi, + IEP.begin(ucell, + this->psi, this->pw_rhod, this->pw_wfc, this->pw_big, @@ -344,7 +348,8 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - IEP.begin(this->psi, + IEP.begin(ucell, + this->psi, this->pw_rhod, this->pw_wfc, this->pw_big, diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index c36775aa07..ef32f041e8 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -56,13 +56,17 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) const std::string cal_type = PARAM.inp.calculation; if (cal_type == "test_memory") { - Cal_Test::test_memory(this->pw_rho, + Cal_Test::test_memory(ucell.nat, + ucell.ntype, + ucell.GGT, + this->pw_rho, this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); } else if (cal_type == "gen_bessel") { Numerical_Descriptor nc; - nc.output_descriptor(this->psi[0], + nc.output_descriptor(ucell, + this->psi[0], PARAM.inp.bessel_descriptor_lmax, PARAM.inp.bessel_descriptor_rcut, PARAM.inp.bessel_descriptor_tolerence, diff --git a/source/module_hamilt_lcao/module_gint/gint_gamma.h b/source/module_hamilt_lcao/module_gint/gint_gamma.h index d27c20c53c..323893cbdf 100644 --- a/source/module_hamilt_lcao/module_gint/gint_gamma.h +++ b/source/module_hamilt_lcao/module_gint/gint_gamma.h @@ -36,7 +36,7 @@ class Gint_Gamma : public Gint //! in gint_gamma_env.cpp //! calcualte the electronic wave functions via grid integral - void cal_env(const double* wfc, double* rho,UnitCell &ucell); + void cal_env(const double* wfc, double* rho,const UnitCell &ucell); //! transfer this->hRGint to Veff::hR void transfer_pvpR(hamilt::HContainer* hR,const UnitCell* ucell); diff --git a/source/module_hamilt_lcao/module_gint/gint_gamma_env.cpp b/source/module_hamilt_lcao/module_gint/gint_gamma_env.cpp index 708fea818a..222c604c73 100644 --- a/source/module_hamilt_lcao/module_gint/gint_gamma_env.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_gamma_env.cpp @@ -6,7 +6,7 @@ #include "module_basis/module_ao/ORB_read.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -void Gint_Gamma::cal_env(const double* wfc, double* rho, UnitCell& ucell) +void Gint_Gamma::cal_env(const double* wfc, double* rho,const UnitCell& ucell) { ModuleBase::TITLE("Grid_Integral", "cal_env"); diff --git a/source/module_hamilt_lcao/module_gint/gint_k.h b/source/module_hamilt_lcao/module_gint/gint_k.h index 21968b3a27..0a609d3ac7 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k.h +++ b/source/module_hamilt_lcao/module_gint/gint_k.h @@ -65,7 +65,7 @@ class Gint_k : public Gint { double* rho, const std::vector>& kvec_c, const std::vector>& kvec_d, - UnitCell& ucell); + const UnitCell& ucell); //------------------------------------------------------ // in gint_k_sparse1.cpp diff --git a/source/module_hamilt_lcao/module_gint/gint_k_env.cpp b/source/module_hamilt_lcao/module_gint/gint_k_env.cpp index b6510a0a51..da212ce4f0 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_env.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_env.cpp @@ -13,7 +13,7 @@ void Gint_k::cal_env_k(int ik, double* rho, const std::vector>& kvec_c, const std::vector>& kvec_d, - UnitCell& ucell) + const UnitCell& ucell) { ModuleBase::TITLE("Gint_k", "cal_env_k"); ModuleBase::timer::tick("Gint_k", "cal_env_k"); diff --git a/source/module_io/berryphase.cpp b/source/module_io/berryphase.cpp index 71b173b095..a8b5eef5b3 100644 --- a/source/module_io/berryphase.cpp +++ b/source/module_io/berryphase.cpp @@ -40,12 +40,12 @@ void berryphase::get_occupation_bands() } #ifdef __LCAO -void berryphase::lcao_init(const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb) +void berryphase::lcao_init(const UnitCell& ucell, const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb) { ModuleBase::TITLE("berryphase", "lcao_init"); - lcao_method.init(grid_tech, kv.get_nkstot(), orb); - lcao_method.cal_R_number(); - lcao_method.cal_orb_overlap(); + lcao_method.init(ucell,grid_tech, kv.get_nkstot(), orb); + lcao_method.cal_R_number(ucell); + lcao_method.cal_orb_overlap(ucell); return; } #endif @@ -232,7 +232,8 @@ void berryphase::set_kpoints(const K_Vectors& kv, const int direction) } #include "../module_base/complexmatrix.h" -double berryphase::stringPhase(int index_str, +double berryphase::stringPhase(const UnitCell& ucell, + int index_str, int nbands, const int npwx, const psi::Psi>* psi_in, @@ -347,7 +348,7 @@ double berryphase::stringPhase(int index_str, if (PARAM.inp.nspin != 4) { // std::complex my_det = lcao_method.det_berryphase(ik_1,ik_2,dk,nbands); - zeta = zeta * lcao_method.det_berryphase(ik_1, ik_2, dk, nbands, *(this->paraV), psi_in, kv); + zeta = zeta * lcao_method.det_berryphase(ucell,ik_1, ik_2, dk, nbands, *(this->paraV), psi_in, kv); // test by jingan // GlobalV::ofs_running << "methon 1: det = " << my_det << std::endl; // test by jingan @@ -392,7 +393,8 @@ double berryphase::stringPhase(int index_str, return log(zeta).imag(); } -void berryphase::Berry_Phase(int nbands, +void berryphase::Berry_Phase(const UnitCell& ucell, + int nbands, double& pdl_elec_tot, int& mod_elec_tot, const int npwx, @@ -421,7 +423,7 @@ void berryphase::Berry_Phase(int nbands, for (int istring = 0; istring < total_string; istring++) { - phik[istring] = stringPhase(istring, nbands, npwx, psi_in, rhopw, wfcpw, kv); + phik[istring] = stringPhase(ucell,istring, nbands, npwx, psi_in, rhopw, wfcpw, kv); // transfer phase to complex number cphik[istring] = std::complex(cos(phik[istring]), sin(phik[istring])); cave = cave + std::complex(wistring[istring], 0.0) * cphik[istring]; @@ -463,7 +465,8 @@ void berryphase::Berry_Phase(int nbands, // GlobalV::ofs_running << "Berry_Phase end " << std::endl; } -void berryphase::Macroscopic_polarization(const int npwx, +void berryphase::Macroscopic_polarization(const UnitCell& ucell, + const int npwx, const psi::Psi>* psi_in, const ModulePW::PW_Basis* rhopw, const ModulePW::PW_Basis_K* wfcpw, @@ -486,30 +489,30 @@ void berryphase::Macroscopic_polarization(const int npwx, double polarization_ion[3]; // means three lattice vector directions R1，R2，R3 ModuleBase::GlobalFunc::ZEROS(polarization_ion, 3); // reciprocal lattice - ModuleBase::Vector3 rcell_1(GlobalC::ucell.G.e11, GlobalC::ucell.G.e12, GlobalC::ucell.G.e13); - ModuleBase::Vector3 rcell_2(GlobalC::ucell.G.e21, GlobalC::ucell.G.e22, GlobalC::ucell.G.e23); - ModuleBase::Vector3 rcell_3(GlobalC::ucell.G.e31, GlobalC::ucell.G.e32, GlobalC::ucell.G.e33); - // int *mod_ion = new int[GlobalC::ucell.nat]; - std::vector mod_ion(GlobalC::ucell.nat); - // double *pdl_ion_R1 = new double[GlobalC::ucell.nat]; - std::vector pdl_ion_R1(GlobalC::ucell.nat); - // double *pdl_ion_R2 = new double[GlobalC::ucell.nat]; - std::vector pdl_ion_R2(GlobalC::ucell.nat); - // double *pdl_ion_R3 = new double[GlobalC::ucell.nat]; - std::vector pdl_ion_R3(GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(mod_ion.data(), GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pdl_ion_R1.data(), GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pdl_ion_R2.data(), GlobalC::ucell.nat); - ModuleBase::GlobalFunc::ZEROS(pdl_ion_R3.data(), GlobalC::ucell.nat); + ModuleBase::Vector3 rcell_1(ucell.G.e11, ucell.G.e12, ucell.G.e13); + ModuleBase::Vector3 rcell_2(ucell.G.e21, ucell.G.e22, ucell.G.e23); + ModuleBase::Vector3 rcell_3(ucell.G.e31, ucell.G.e32, ucell.G.e33); + // int *mod_ion = new int[ucell.nat]; + std::vector mod_ion(ucell.nat); + // double *pdl_ion_R1 = new double[ucell.nat]; + std::vector pdl_ion_R1(ucell.nat); + // double *pdl_ion_R2 = new double[ucell.nat]; + std::vector pdl_ion_R2(ucell.nat); + // double *pdl_ion_R3 = new double[ucell.nat]; + std::vector pdl_ion_R3(ucell.nat); + ModuleBase::GlobalFunc::ZEROS(mod_ion.data(), ucell.nat); + ModuleBase::GlobalFunc::ZEROS(pdl_ion_R1.data(), ucell.nat); + ModuleBase::GlobalFunc::ZEROS(pdl_ion_R2.data(), ucell.nat); + ModuleBase::GlobalFunc::ZEROS(pdl_ion_R3.data(), ucell.nat); bool lodd = false; int atom_index = 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++) { // should consider fractional electron number - if (int(GlobalC::ucell.atoms[it].ncpp.zv) % 2 == 1) + if (int(ucell.atoms[it].ncpp.zv) % 2 == 1) { mod_ion[atom_index] = 1; lodd = true; @@ -524,18 +527,18 @@ void berryphase::Macroscopic_polarization(const int npwx, } atom_index = 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++) { - pdl_ion_R1[atom_index] = GlobalC::ucell.atoms[it].ncpp.zv * (GlobalC::ucell.atoms[it].tau[ia] * rcell_1); - pdl_ion_R2[atom_index] = GlobalC::ucell.atoms[it].ncpp.zv * (GlobalC::ucell.atoms[it].tau[ia] * rcell_2); - pdl_ion_R3[atom_index] = GlobalC::ucell.atoms[it].ncpp.zv * (GlobalC::ucell.atoms[it].tau[ia] * rcell_3); + pdl_ion_R1[atom_index] = ucell.atoms[it].ncpp.zv * (ucell.atoms[it].tau[ia] * rcell_1); + pdl_ion_R2[atom_index] = ucell.atoms[it].ncpp.zv * (ucell.atoms[it].tau[ia] * rcell_2); + pdl_ion_R3[atom_index] = ucell.atoms[it].ncpp.zv * (ucell.atoms[it].tau[ia] * rcell_3); atom_index++; } } - for (int i = 0; i < GlobalC::ucell.nat; i++) + for (int i = 0; i < ucell.nat; i++) { if (mod_ion[i] == 1) { @@ -594,12 +597,12 @@ void berryphase::Macroscopic_polarization(const int npwx, set_kpoints(kv, direction); double pdl_elec_tot = 0.0; int mod_elec_tot = 0; - Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); + Berry_Phase(ucell,occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); - const double rmod = GlobalC::ucell.a1.norm() * GlobalC::ucell.lat0; + const double rmod = ucell.a1.norm() * ucell.lat0; const double unit1 = rmod; - const double unit2 = rmod / GlobalC::ucell.omega; - const double unit3 = (rmod / GlobalC::ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); + const double unit2 = rmod / ucell.omega; + const double unit3 = (rmod / ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); GlobalV::ofs_running << " VALUES OF POLARIZATION" << std::endl; GlobalV::ofs_running << std::endl; @@ -613,7 +616,7 @@ void berryphase::Macroscopic_polarization(const int npwx, // calculate total polarization,add electron part and ions part double total_polarization = pdl_elec_tot + polarization_ion[0]; - ModuleBase::Vector3 polarization_xyz = GlobalC::ucell.a1; + ModuleBase::Vector3 polarization_xyz = ucell.a1; polarization_xyz.normalize(); polarization_xyz = total_polarization * polarization_xyz; @@ -640,12 +643,12 @@ void berryphase::Macroscopic_polarization(const int npwx, set_kpoints(kv, direction); double pdl_elec_tot = 0.0; int mod_elec_tot = 0; - Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); + Berry_Phase(ucell,occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); - const double rmod = GlobalC::ucell.a2.norm() * GlobalC::ucell.lat0; + const double rmod = ucell.a2.norm() * ucell.lat0; const double unit1 = rmod; - const double unit2 = rmod / GlobalC::ucell.omega; - const double unit3 = (rmod / GlobalC::ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); + const double unit2 = rmod / ucell.omega; + const double unit3 = (rmod / ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); GlobalV::ofs_running << " VALUES OF POLARIZATION" << std::endl; GlobalV::ofs_running << std::endl; @@ -659,7 +662,7 @@ void berryphase::Macroscopic_polarization(const int npwx, // calculate total polarization,add electron part and ions part double total_polarization = pdl_elec_tot + polarization_ion[1]; - ModuleBase::Vector3 polarization_xyz = GlobalC::ucell.a2; + ModuleBase::Vector3 polarization_xyz = ucell.a2; polarization_xyz.normalize(); polarization_xyz = total_polarization * polarization_xyz; @@ -686,12 +689,12 @@ void berryphase::Macroscopic_polarization(const int npwx, set_kpoints(kv, direction); double pdl_elec_tot = 0.0; int mod_elec_tot = 0; - Berry_Phase(occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); + Berry_Phase(ucell,occ_nbands, pdl_elec_tot, mod_elec_tot, npwx, psi_in, rhopw, wfcpw, kv); - const double rmod = GlobalC::ucell.a3.norm() * GlobalC::ucell.lat0; + const double rmod = ucell.a3.norm() * ucell.lat0; const double unit1 = rmod; - const double unit2 = rmod / GlobalC::ucell.omega; - const double unit3 = (rmod / GlobalC::ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); + const double unit2 = rmod / ucell.omega; + const double unit3 = (rmod / ucell.omega) * (1.60097e-19 / pow(5.29177e-11, 2)); GlobalV::ofs_running << " VALUES OF POLARIZATION" << std::endl; GlobalV::ofs_running << std::endl; @@ -705,7 +708,7 @@ void berryphase::Macroscopic_polarization(const int npwx, // calculate total polarization,add electron part and ions part double total_polarization = pdl_elec_tot + polarization_ion[2]; - ModuleBase::Vector3 polarization_xyz = GlobalC::ucell.a3; + ModuleBase::Vector3 polarization_xyz = ucell.a3; polarization_xyz.normalize(); polarization_xyz = total_polarization * polarization_xyz; diff --git a/source/module_io/berryphase.h b/source/module_io/berryphase.h index 0e226b7cad..7a3707c711 100644 --- a/source/module_io/berryphase.h +++ b/source/module_io/berryphase.h @@ -37,11 +37,12 @@ class berryphase void get_occupation_bands(); #ifdef __LCAO - void lcao_init(const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb); + void lcao_init(const UnitCell& ucell, const K_Vectors& kv, const Grid_Technique& grid_tech, const LCAO_Orbitals& orb); #endif void set_kpoints(const K_Vectors& kv, const int direction); - double stringPhase(int index_str, + double stringPhase(const UnitCell& ucell, + int index_str, int nbands, const int npwx, const psi::Psi>* psi_in, @@ -49,7 +50,8 @@ class berryphase const ModulePW::PW_Basis_K* wfcpw, const K_Vectors& kv); - void Berry_Phase(int nbands, + void Berry_Phase(const UnitCell& ucell, + int nbands, double& pdl_elec_tot, int& mod_elec_tot, const int npwx, @@ -58,7 +60,8 @@ class berryphase const ModulePW::PW_Basis_K* wfcpw, const K_Vectors& kv); - void Macroscopic_polarization(const int npwx, + void Macroscopic_polarization(const UnitCell& ucell, + const int npwx, const psi::Psi* psi_in, const ModulePW::PW_Basis* rhopw, const ModulePW::PW_Basis_K* wfcpw, @@ -66,7 +69,8 @@ class berryphase { throw std::logic_error("berry phase supports only multi-k"); }; - void Macroscopic_polarization(const int npwx, + void Macroscopic_polarization(const UnitCell& ucell, + const int npwx, const psi::Psi>* psi_in, const ModulePW::PW_Basis* rhopw, const ModulePW::PW_Basis_K* wfcpw, diff --git a/source/module_io/cal_r_overlap_R.cpp b/source/module_io/cal_r_overlap_R.cpp index 0d1da466b8..b43c181576 100644 --- a/source/module_io/cal_r_overlap_R.cpp +++ b/source/module_io/cal_r_overlap_R.cpp @@ -14,7 +14,8 @@ cal_r_overlap_R::~cal_r_overlap_R() { } -void cal_r_overlap_R::initialize_orb_table(const LCAO_Orbitals& orb) +void cal_r_overlap_R::initialize_orb_table(const UnitCell& ucell, + const LCAO_Orbitals& orb) { int Lmax_used = 0; int Lmax = 0; @@ -28,7 +29,7 @@ void cal_r_overlap_R::initialize_orb_table(const LCAO_Orbitals& orb) for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -53,7 +54,8 @@ void cal_r_overlap_R::initialize_orb_table(const LCAO_Orbitals& orb) MGT.init_Gaunt(Lmax); } -void cal_r_overlap_R::construct_orbs_and_orb_r(const LCAO_Orbitals& orb) +void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, + const LCAO_Orbitals& orb) { int orb_r_ntype = 0; int mat_Nr = orb.Phi[0].PhiLN(0, 0).getNr(); @@ -201,13 +203,13 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const LCAO_Orbitals& orb) iw2im.resize(PARAM.globalv.nlocal); int iw = 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++) { - for (int iL = 0; iL < GlobalC::ucell.atoms[it].nwl + 1; iL++) + for (int iL = 0; iL < ucell.atoms[it].nwl + 1; iL++) { - for (int iN = 0; iN < GlobalC::ucell.atoms[it].l_nchi[iL]; iN++) + for (int iN = 0; iN < ucell.atoms[it].l_nchi[iL]; iN++) { for (int im = 0; im < (2 * iL + 1); im++) { @@ -224,20 +226,20 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const LCAO_Orbitals& orb) } } -void cal_r_overlap_R::init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) +void cal_r_overlap_R::init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) { ModuleBase::TITLE("cal_r_overlap_R", "init"); ModuleBase::timer::tick("cal_r_overlap_R", "init"); this->ParaV = &pv; - initialize_orb_table(orb); - construct_orbs_and_orb_r(orb); + initialize_orb_table(ucell,orb); + construct_orbs_and_orb_r(ucell,orb); ModuleBase::timer::tick("cal_r_overlap_R", "init"); return; } -void cal_r_overlap_R::out_rR(const int& istep) +void cal_r_overlap_R::out_rR(const UnitCell& ucell, const int& istep) { ModuleBase::TITLE("cal_r_overlap_R", "out_rR"); ModuleBase::timer::tick("cal_r_overlap_R", "out_rR"); @@ -296,7 +298,7 @@ void cal_r_overlap_R::out_rR(const int& istep) int dRy = R_coor.y; int dRz = R_coor.z; - ModuleBase::Vector3 R_car = ModuleBase::Vector3(dRx, dRy, dRz) * GlobalC::ucell.latvec; + ModuleBase::Vector3 R_car = ModuleBase::Vector3(dRx, dRy, dRz) * ucell.latvec; int ir, ic; for (int iw1 = 0; iw1 < PARAM.globalv.nlocal; iw1++) @@ -331,8 +333,8 @@ void cal_r_overlap_R::out_rR(const int& istep) int im2 = iw2im[orb_index_col]; ModuleBase::Vector3 r_distance - = (GlobalC::ucell.atoms[it2].tau[ia2] - GlobalC::ucell.atoms[it1].tau[ia1] + R_car) - * GlobalC::ucell.lat0; + = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) + * ucell.lat0; double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, r_distance, @@ -365,7 +367,7 @@ void cal_r_overlap_R::out_rR(const int& istep) ModuleBase::Vector3 temp_prp = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) - + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o; + + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; if (std::abs(temp_prp.x) > sparse_threshold) { @@ -506,7 +508,7 @@ void cal_r_overlap_R::out_rR(const int& istep) return; } -void cal_r_overlap_R::out_rR_other(const int& istep, const std::set>& output_R_coor) +void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, const int& istep, const std::set>& output_R_coor) { ModuleBase::TITLE("cal_r_overlap_R", "out_rR_other"); ModuleBase::timer::tick("cal_r_overlap_R", "out_rR_other"); @@ -571,7 +573,7 @@ void cal_r_overlap_R::out_rR_other(const int& istep, const std::set R_car = ModuleBase::Vector3(dRx, dRy, dRz) * GlobalC::ucell.latvec; + ModuleBase::Vector3 R_car = ModuleBase::Vector3(dRx, dRy, dRz) * ucell.latvec; int ir = 0; int ic = 0; @@ -607,8 +609,8 @@ void cal_r_overlap_R::out_rR_other(const int& istep, const std::set r_distance - = (GlobalC::ucell.atoms[it2].tau[ia2] - GlobalC::ucell.atoms[it1].tau[ia1] + R_car) - * GlobalC::ucell.lat0; + = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) + * ucell.lat0; double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, r_distance, @@ -641,7 +643,7 @@ void cal_r_overlap_R::out_rR_other(const int& istep, const std::set temp_prp = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) - + GlobalC::ucell.atoms[it1].tau[ia1] * GlobalC::ucell.lat0 * overlap_o; + + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; if (std::abs(temp_prp.x) > sparse_threshold) { diff --git a/source/module_io/cal_r_overlap_R.h b/source/module_io/cal_r_overlap_R.h index 0b2ff6b6fd..ae23718148 100644 --- a/source/module_io/cal_r_overlap_R.h +++ b/source/module_io/cal_r_overlap_R.h @@ -13,7 +13,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb21.h" #include "module_hamilt_lcao/hamilt_lcaodft/center2_orb.h" #include "single_R_io.h" - +#include "module_cell/unitcell.h" #include #include #include @@ -30,13 +30,13 @@ class cal_r_overlap_R double sparse_threshold = 1e-10; bool binary = false; - void init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); - void out_rR(const int& istep); - void out_rR_other(const int& istep, const std::set>& output_R_coor); + void init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); + void out_rR(const UnitCell& ucell, const int& istep); + void out_rR_other(const UnitCell& ucell, const int& istep, const std::set>& output_R_coor); private: - void initialize_orb_table(const LCAO_Orbitals& orb); - void construct_orbs_and_orb_r(const LCAO_Orbitals& orb); + void initialize_orb_table(const UnitCell& ucell, const LCAO_Orbitals& orb); + void construct_orbs_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb); std::vector iw2ia; std::vector iw2iL; diff --git a/source/module_io/cal_test.cpp b/source/module_io/cal_test.cpp index 680ba6028f..b41c38a2ba 100644 --- a/source/module_io/cal_test.cpp +++ b/source/module_io/cal_test.cpp @@ -48,16 +48,18 @@ double Cal_Test::meigts123=0.0; double Cal_Test::mtot=0.0; -void Cal_Test::test_memory( - const ModulePW::PW_Basis* rhopw, - const ModulePW::PW_Basis_K* wfcpw, - const std::string chr_mixing_mode, - const int chr_mixing_ndim) +void Cal_Test::test_memory(const int nat, + const int ntype, + const ModuleBase::Matrix3& GGT, + const ModulePW::PW_Basis* rhopw, + const ModulePW::PW_Basis_K* wfcpw, + const std::string chr_mixing_mode, + const int chr_mixing_ndim) { ModuleBase::TITLE("Cal_Test","test_memory"); - const int ngmw = Cal_Test::cal_np(wfcpw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz); - const int ngmc = Cal_Test::cal_np(rhopw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz); + const int ngmw = Cal_Test::cal_np(GGT,wfcpw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz); + const int ngmc = Cal_Test::cal_np(GGT,rhopw->ggecut, rhopw->nx, rhopw->ny, rhopw->nz); // const int ecut_wfc = INPUT.ecutwfc; // const int ecut_chg = INPUT.ecutrho; @@ -65,7 +67,7 @@ void Cal_Test::test_memory( // const int ngmw = Cal_Test::cal_np(ecut_wfc, rhopw->nx, rhopw->ny, rhopw->nz); // const int ngmc = Cal_Test::cal_np(ecut_chg, rhopw->nx, rhopw->ny, rhopw->nz); - std::cout << " number of atoms = " << GlobalC::ucell.nat << std::endl; + std::cout << " number of atoms = " << nat << std::endl; std::cout << " plane wave number for wave functions = " << ngmw << std::endl; std::cout << " plane wave number for chage density = " << ngmc << std::endl; @@ -108,8 +110,8 @@ void Cal_Test::test_memory( mig2fftc = ModuleBase::Memory::calculate_mem( ngmc , "int"); mgg = ModuleBase::Memory::calculate_mem( ngmc, "double"); mig123 = ModuleBase::Memory::calculate_mem( ngmc*3, "int"); - mstrucFac = ModuleBase::Memory::calculate_mem( GlobalC::ucell.ntype*ngmc, "cdouble"); - meigts123 = ModuleBase::Memory::calculate_mem( GlobalC::ucell.nat * (2*rhopw->nx+1+2*rhopw->ny+1+2*rhopw->nz+1), "cdouble"); + mstrucFac = ModuleBase::Memory::calculate_mem( ntype*ngmc, "cdouble"); + meigts123 = ModuleBase::Memory::calculate_mem( nat * (2*rhopw->nx+1+2*rhopw->ny+1+2*rhopw->nz+1), "cdouble"); //(3) Memory for H,S matrix. std::cout << " NLOCAL = " << PARAM.globalv.nlocal << std::endl; @@ -130,7 +132,7 @@ void Cal_Test::test_memory( // print_mem(8); // print_mem(16); -// if(GlobalC::ucell.nat > 200) +// if(nat > 200) // { // print_mem(32); // print_mem(64); @@ -140,7 +142,11 @@ void Cal_Test::test_memory( } //! compute the number of plane waves -int Cal_Test::cal_np(const double &ggcut, const int &n1, const int &n2, const int &n3) +int Cal_Test::cal_np(const ModuleBase::Matrix3& GGT, + const double &ggcut, + const int &n1, + const int &n2, + const int &n3) { /* @@ -170,7 +176,7 @@ int Cal_Test::cal_np(const double &ggcut, const int &n1, const int &n2, const in { ModuleBase::Vector3 f(i,j,k); // g2= |f|^2 in the unit of (2Pi/lat0)^2 - double g2 = f * (GlobalC::ucell.GGT * f); + double g2 = f * (GGT * f); // gcut is from input. if (g2 <= ggcut) diff --git a/source/module_io/cal_test.h b/source/module_io/cal_test.h index be272e914b..1134ec06ae 100644 --- a/source/module_io/cal_test.h +++ b/source/module_io/cal_test.h @@ -5,12 +5,19 @@ namespace Cal_Test { - void test_memory(const ModulePW::PW_Basis* rhopw, - const ModulePW::PW_Basis_K* wfcpw, - const std::string chr_mixing_mode, - const int chr_mixing_ndim); + void test_memory(const int nat, + const int ntype, + const ModuleBase::Matrix3& GGT, + const ModulePW::PW_Basis* rhopw, + const ModulePW::PW_Basis_K* wfcpw, + const std::string chr_mixing_mode, + const int chr_mixing_ndim); - int cal_np(const double &ggcut, const int &n1, const int &n2, const int &n3); + int cal_np(const ModuleBase::Matrix3& GGT, + const double &ggcut, + const int &n1, + const int &n2, + const int &n3); void print_mem(const int &nproc); diff --git a/source/module_io/dipole_io.h b/source/module_io/dipole_io.h index d9fbcfc74f..16c33e6fed 100644 --- a/source/module_io/dipole_io.h +++ b/source/module_io/dipole_io.h @@ -7,7 +7,8 @@ namespace ModuleIO { -void write_dipole(const double* rho_save, +void write_dipole(const UnitCell& ucell, + const double* rho_save, const ModulePW::PW_Basis* rhopw, const int& is, const int& istep, diff --git a/source/module_io/dos_nao.cpp b/source/module_io/dos_nao.cpp index fb993cc201..894673bfe1 100644 --- a/source/module_io/dos_nao.cpp +++ b/source/module_io/dos_nao.cpp @@ -6,18 +6,18 @@ namespace ModuleIO { /// @brief manege the output of dos in numerical atomic basis case -/// @param[in] psi -/// @param[in] pv -/// @param[in] ekb -/// @param[in] wg -/// @param[in] dos_edelta_ev -/// @param[in] dos_scale -/// @param[in] dos_sigma -/// @param[in] kv -/// @param[in] Pkpoints -/// @param[in] ucell -/// @param[in] eferm -/// @param[in] nbands + /// @param[in] psi + /// @param[in] pv + /// @param[in] ekb + /// @param[in] wg + /// @param[in] dos_edelta_ev + /// @param[in] dos_scale + /// @param[in] dos_sigma + /// @param[in] kv + /// @param[in] Pkpoints + /// @param[in] ucell + /// @param[in] eferm + /// @param[in] nbands template void out_dos_nao( const psi::Psi* psi, @@ -36,7 +36,7 @@ namespace ModuleIO { ModuleBase::TITLE("Module_IO", "out_dos_nao"); - write_dos_lcao(psi, pv, ekb, wg, dos_edelta_ev, dos_scale, dos_sigma, kv, p_ham); + write_dos_lcao(ucell,psi, pv, ekb, wg, dos_edelta_ev, dos_scale, dos_sigma, kv, p_ham); int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1; if (PARAM.inp.out_dos == 3) diff --git a/source/module_io/get_wf_lcao.cpp b/source/module_io/get_wf_lcao.cpp index d1257de1bb..7197bb8292 100644 --- a/source/module_io/get_wf_lcao.cpp +++ b/source/module_io/get_wf_lcao.cpp @@ -19,7 +19,8 @@ IState_Envelope::~IState_Envelope() } // For gamma_only -void IState_Envelope::begin(const psi::Psi* psid, +void IState_Envelope::begin(const UnitCell& ucell, + const psi::Psi* psid, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, @@ -113,7 +114,7 @@ void IState_Envelope::begin(const psi::Psi* psid, } #endif - gg.cal_env(wfc_gamma_grid[is][ib], pes_->charge->rho[is], GlobalC::ucell); + gg.cal_env(wfc_gamma_grid[is][ib], pes_->charge->rho[is], ucell); pes_->charge->save_rho_before_sum_band(); @@ -128,7 +129,7 @@ void IState_Envelope::begin(const psi::Psi* psid, 0, ss.str(), ef_tmp, - &(GlobalC::ucell)); + &(ucell)); } } } @@ -176,7 +177,7 @@ void IState_Envelope::begin(const psi::Psi* psid, } #endif - gg.cal_env(wfc_gamma_grid[is][ib], pes_->charge->rho[is], GlobalC::ucell); + gg.cal_env(wfc_gamma_grid[is][ib], pes_->charge->rho[is], ucell); pes_->charge->save_rho_before_sum_band(); @@ -210,7 +211,7 @@ void IState_Envelope::begin(const psi::Psi* psid, 0, ss_real.str(), ef_tmp, - &(GlobalC::ucell)); + &(ucell)); // Output imaginary part std::stringstream ss_imag; @@ -222,7 +223,7 @@ void IState_Envelope::begin(const psi::Psi* psid, 0, ss_imag.str(), ef_tmp, - &(GlobalC::ucell)); + &(ucell)); } } } @@ -238,7 +239,7 @@ void IState_Envelope::begin(const psi::Psi* psid, if (out_wfc_r) { - ModuleIO::write_psi_r_1(psi_g, pw_wfc, "wfc_realspace", false, kv); + ModuleIO::write_psi_r_1(ucell,psi_g, pw_wfc, "wfc_realspace", false, kv); } for (int is = 0; is < nspin; ++is) @@ -253,7 +254,8 @@ void IState_Envelope::begin(const psi::Psi* psid, } // For multi-k -void IState_Envelope::begin(const psi::Psi>* psi, +void IState_Envelope::begin(const UnitCell& ucell, + const psi::Psi>* psi, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, @@ -348,7 +350,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, } #endif // deal with NSPIN=4 - gk.cal_env_k(ik, wfc_k_grid[ik][ib], pes_->charge->rho[ispin], kv.kvec_c, kv.kvec_d, GlobalC::ucell); + gk.cal_env_k(ik, wfc_k_grid[ik][ib], pes_->charge->rho[ispin], kv.kvec_c, kv.kvec_d, ucell); std::stringstream ss; ss << global_out_dir << "BAND" << ib + 1 << "_k_" << ik + 1 << "_s_" << ispin + 1 << "_ENV.cube"; @@ -361,7 +363,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, 0, ss.str(), ef_tmp, - &(GlobalC::ucell), + &(ucell), 3, 1); @@ -386,7 +388,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, } if (out_wf_r) { - ModuleIO::write_psi_r_1(psi_g, pw_wfc, "wfc_realspace", false, kv); + ModuleIO::write_psi_r_1(ucell,psi_g, pw_wfc, "wfc_realspace", false, kv); } std::cout << " Outputting real-space wave functions in cube format..." << std::endl; @@ -429,7 +431,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, 0, ss_real.str(), ef_tmp, - &(GlobalC::ucell)); + &(ucell)); // Output imaginary part std::stringstream ss_imag; @@ -442,7 +444,7 @@ void IState_Envelope::begin(const psi::Psi>* psi, 0, ss_imag.str(), ef_tmp, - &(GlobalC::ucell)); + &(ucell)); } } } diff --git a/source/module_io/get_wf_lcao.h b/source/module_io/get_wf_lcao.h index b4c8c11d56..df5ed1e735 100644 --- a/source/module_io/get_wf_lcao.h +++ b/source/module_io/get_wf_lcao.h @@ -17,7 +17,8 @@ class IState_Envelope ~IState_Envelope(); /// For gamma_only - void begin(const psi::Psi* psid, + void begin(const UnitCell& ucell, + const psi::Psi* psid, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, @@ -36,7 +37,8 @@ class IState_Envelope const std::string& global_out_dir); /// tmp, delete after Gint is refactored. - void begin(const psi::Psi* psid, + void begin(const UnitCell& ucell, + const psi::Psi* psid, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, @@ -58,7 +60,8 @@ class IState_Envelope }; /// For multi-k - void begin(const psi::Psi>* psi, + void begin(const UnitCell& ucell, + const psi::Psi>* psi, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, @@ -77,7 +80,8 @@ class IState_Envelope const std::string& global_out_dir); /// tmp, delete after Gint is refactored. - void begin(const psi::Psi>* psi, + void begin(const UnitCell& ucell, + const psi::Psi>* psi, const ModulePW::PW_Basis* pw_rhod, const ModulePW::PW_Basis_K* pw_wfc, const ModulePW::PW_Basis_Big* pw_big, diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index ac93d06bd9..5b807620c8 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -255,7 +255,7 @@ void Input_Conv::Convert() GlobalC::dftu.U0 = PARAM.globalv.hubbard_u; if (PARAM.globalv.uramping > 0.01) { - ModuleBase::GlobalFunc::ZEROS(GlobalC::dftu.U.data(), GlobalC::ucell.ntype); + ModuleBase::GlobalFunc::ZEROS(GlobalC::dftu.U.data(), PARAM.inp.ntype); } } #endif diff --git a/source/module_io/numerical_descriptor.cpp b/source/module_io/numerical_descriptor.cpp index aea1840e11..aa879d40af 100644 --- a/source/module_io/numerical_descriptor.cpp +++ b/source/module_io/numerical_descriptor.cpp @@ -25,7 +25,7 @@ Numerical_Descriptor::~Numerical_Descriptor() } -void Numerical_Descriptor::output_descriptor(const psi::Psi> &psi, const int &lmax_in, const double &rcut_in, const double &tol_in, const int nks_in) +void Numerical_Descriptor::output_descriptor(const UnitCell& ucell, const psi::Psi> &psi, const int &lmax_in, const double &rcut_in, const double &tol_in, const int nks_in) { ModuleBase::TITLE("Numerical_Descriptor","output_descriptor"); ModuleBase::GlobalFunc::NEW_PART("DeepKS descriptor: D_{Inl}"); @@ -49,16 +49,16 @@ void Numerical_Descriptor::output_descriptor(const psi::Psi this->bessel_basis.init( false, std::stod(PARAM.inp.bessel_descriptor_ecut), - GlobalC::ucell.ntype, + ucell.ntype, this->lmax, PARAM.inp.bessel_descriptor_smooth, PARAM.inp.bessel_descriptor_sigma, rcut_in, tol_in, - GlobalC::ucell + ucell ); this->nmax = Numerical_Descriptor::bessel_basis.get_ecut_number(); - this->init_mu_index(); + this->init_mu_index(ucell); this->init_label = true; assert(nmax>0); @@ -129,10 +129,10 @@ void Numerical_Descriptor::output_descriptor(const psi::Psi // 5. Generate descriptors for each atom //------------------------------------- - for (int it=0; it // for each 'lmax' we have 'n' up to 'ecut_number' this->generate_descriptor(overlap_Q1, overlap_Q2, it ,ia, d, nd); - ofs << GlobalC::ucell.atoms[it].label << " atom_index " << ia+1 << " n_descriptor " << nd << std::endl; + ofs << ucell.atoms[it].label << " atom_index " << ia+1 << " n_descriptor " << nd << std::endl; for(int id=0; id0 && id%8==0) ofs << std::endl; @@ -264,7 +264,7 @@ void Numerical_Descriptor::jlq3d_overlap( GlobalV::ofs_running << " OUTPUT THE OVERLAP BETWEEN SPHERICAL BESSEL FUNCTIONS AND BLOCH WAVE FUNCTIONS" << std::endl; GlobalV::ofs_running << " Q = < J_it_ia_il_in_im | Psi_n, k > " << std::endl; - const double normalization = (4 * ModuleBase::PI) / sqrt(GlobalC::ucell.omega);// Peize Lin add normalization + const double normalization = (4 * ModuleBase::PI) / sqrt(ucell.omega);// Peize Lin add normalization 2015-12-29 const int total_lm = ( this->lmax + 1) * ( this->lmax + 1); @@ -286,15 +286,15 @@ std::endl; GlobalV::ofs_running << " Q = < J_it_ia_il_in_im | Psi_n, k > " << st double *flq = new double[np]; std::complex overlapQ = ModuleBase::ZERO; - for (int T1 = 0; T1 < GlobalC::ucell.ntype; T1++) + for (int T1 = 0; T1 < ucell.ntype; T1++) { - for (int I1 = 0; I1 < GlobalC::ucell.atoms[T1].na; I1++) + for (int I1 = 0; I1 < ucell.atoms[T1].na; I1++) { std::complex *sk = sf.get_sk(ik, T1, I1,wfcpw); for (int L=0; L< lmax+1; L++) { GlobalV::ofs_running << " " << std::setw(5) << ik+1 - << std::setw(8) << GlobalC::ucell.atoms[T1].label + << std::setw(8) << ucell.atoms[T1].label << std::setw(8) << I1+1 << std::setw(8) << L << std::endl; @@ -305,7 +305,7 @@ normalization 2015-12-29 for (int ie=0; ie < nmax; ie++) for (int ig=0; ig> &psi, const int &lmax_in, const double &rcut_in, const double &tol_in, const int nks); // mohan added 2021-01-03 + void output_descriptor(const UnitCell& ucell, const psi::Psi> &psi, const int &lmax_in, const double &rcut_in, const double &tol_in, const int nks); // mohan added 2021-01-03 private: @@ -33,7 +33,7 @@ class Numerical_Descriptor Bessel_Basis bessel_basis; ModuleBase::IntArray *mu_index; - void init_mu_index(void);//mohan added 2021-01-03 + void init_mu_index(const UnitCell& ucell);//mohan added 2021-01-03 // void jlq3d_overlap(ModuleBase::realArray &overlap_Q1, ModuleBase::realArray &overlap_Q2, // const int &ik_ibz, const int &ik, const int &np, const psi::Psi> &psi); diff --git a/source/module_io/output_mat_sparse.cpp b/source/module_io/output_mat_sparse.cpp index b5139876a5..11cd2caf81 100644 --- a/source/module_io/output_mat_sparse.cpp +++ b/source/module_io/output_mat_sparse.cpp @@ -73,14 +73,14 @@ void output_mat_sparse(const bool& out_mat_hsR, if (out_mat_r) { cal_r_overlap_R r_matrix; - r_matrix.init(pv, orb); + r_matrix.init(ucell, pv, orb); if (out_mat_hsR) { - r_matrix.out_rR_other(istep, HS_Arrays.output_R_coor); + r_matrix.out_rR_other(ucell,istep, HS_Arrays.output_R_coor); } else { - r_matrix.out_rR(istep); + r_matrix.out_rR(ucell,istep); } } diff --git a/source/module_io/read_wfc_to_rho.cpp b/source/module_io/read_wfc_to_rho.cpp index bef1a87462..1a9d6c1237 100644 --- a/source/module_io/read_wfc_to_rho.cpp +++ b/source/module_io/read_wfc_to_rho.cpp @@ -134,7 +134,7 @@ void ModuleIO::read_wfc_to_rho(const ModulePW::PW_Basis_K* pw_wfc, Symmetry_rho srho; for (int is = 0; is < nspin; is++) { - srho.begin(is, chg, chg.rhopw, GlobalC::ucell.symm); + srho.begin(is, chg, chg.rhopw, symm); } ModuleBase::timer::tick("ModuleIO", "read_wfc_pw_to_rho"); diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index 31042e73d2..96272a70e7 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -117,7 +117,8 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix, double> ModuleBase::timer::tick("ModuleIO", "cal_tmp_DM"); } -void ModuleIO::write_current(const int istep, +void ModuleIO::write_current(const UnitCell& ucell, + const int istep, const psi::Psi>* psi, const elecstate::ElecState* pelec, const K_Vectors& kv, @@ -133,7 +134,7 @@ void ModuleIO::write_current(const int istep, std::vector>*> current_term = {nullptr, nullptr, nullptr}; if (!TD_Velocity::tddft_velocity) { - cal_current = new TD_current(&GlobalC::ucell, &GlobalC::GridD, pv, orb, intor); + cal_current = new TD_current(&ucell, &GlobalC::GridD, pv, orb, intor); cal_current->calculate_vcomm_r(); cal_current->calculate_grad_term(); for (int dir = 0; dir < 3; dir++) @@ -163,8 +164,8 @@ void ModuleIO::write_current(const int istep, elecstate::cal_dm_psi(DM_real.get_paraV_pointer(), pelec->wg, psi[0], DM_real); // init DMR - DM_real.init_DMR(ra, &GlobalC::ucell); - DM_imag.init_DMR(ra, &GlobalC::ucell); + DM_real.init_DMR(ra, &ucell); + DM_imag.init_DMR(ra, &ucell); int nks = DM_real.get_DMK_nks(); if (PARAM.inp.nspin == 2) @@ -197,27 +198,27 @@ void ModuleIO::write_current(const int istep, #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for (int iat = 0; iat < GlobalC::ucell.nat; iat++) + for (int iat = 0; iat < ucell.nat; iat++) { - const int T1 = GlobalC::ucell.iat2it[iat]; - Atom* atom1 = &GlobalC::ucell.atoms[T1]; - const int I1 = GlobalC::ucell.iat2ia[iat]; + const int T1 = ucell.iat2it[iat]; + Atom* atom1 = &ucell.atoms[T1]; + const int I1 = ucell.iat2ia[iat]; // get iat1 - int iat1 = GlobalC::ucell.itia2iat(T1, I1); + int iat1 = ucell.itia2iat(T1, I1); int irr = pv->nlocstart[iat]; - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + const int start1 = ucell.itiaiw2iwt(T1, I1, 0); for (int cb = 0; cb < ra.na_each[iat]; ++cb) { const int T2 = ra.info[iat][cb][3]; const int I2 = ra.info[iat][cb][4]; - const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + const int start2 = ucell.itiaiw2iwt(T2, I2, 0); - Atom* atom2 = &GlobalC::ucell.atoms[T2]; + Atom* atom2 = &ucell.atoms[T2]; // get iat2 - int iat2 = GlobalC::ucell.itia2iat(T2, I2); + int iat2 = ucell.itia2iat(T2, I2); double Rx = ra.info[iat][cb][0]; double Ry = ra.info[iat][cb][1]; double Rz = ra.info[iat][cb][2]; diff --git a/source/module_io/td_current_io.h b/source/module_io/td_current_io.h index 74acc03134..40aed95214 100644 --- a/source/module_io/td_current_io.h +++ b/source/module_io/td_current_io.h @@ -10,7 +10,8 @@ namespace ModuleIO { #ifdef __LCAO /// @brief func to output current, only used in tddft -void write_current(const int istep, +void write_current(const UnitCell& ucell, + const int istep, const psi::Psi>* psi, const elecstate::ElecState* pelec, const K_Vectors& kv, diff --git a/source/module_io/test/bessel_basis_test.cpp b/source/module_io/test/bessel_basis_test.cpp index e681766090..dc20a25d96 100644 --- a/source/module_io/test/bessel_basis_test.cpp +++ b/source/module_io/test/bessel_basis_test.cpp @@ -514,8 +514,8 @@ TEST_F(TestBesselBasis, PolynomialInterpolationTest) { /* manipulate Bessel_Basis::init_Faln function because for(int it=0; it>* psi = nullptr; Charge chg; - + ModuleSymmetry::Symmetry symm; virtual void SetUp() { wfcpw = new ModulePW::PW_Basis_K; @@ -231,7 +230,7 @@ TEST_F(ReadWfcRhoTest, ReadWfcRho) ModuleIO::write_wfc_pw("WAVEFUNC", *psi, *kv, wfcpw); // Read the wave functions to charge density - ModuleIO::read_wfc_to_rho(wfcpw, GlobalC::ucell.symm, nkstot, kv->isk, chg); + ModuleIO::read_wfc_to_rho(wfcpw, symm, nkstot, kv->isk, chg); // compare the charge density for (int ir = 0; ir < rhopw->nrxx; ++ir) diff --git a/source/module_io/to_wannier90.cpp b/source/module_io/to_wannier90.cpp index aa54f6738d..49a068637e 100644 --- a/source/module_io/to_wannier90.cpp +++ b/source/module_io/to_wannier90.cpp @@ -70,7 +70,7 @@ void toWannier90::calculate() { } -void toWannier90::read_nnkp(const K_Vectors& kv) +void toWannier90::read_nnkp(const UnitCell& ucell, const K_Vectors& kv) { // read *.nnkp file GlobalV::ofs_running << "reading the " << wannier_file_name << ".nnkp file." << std::endl; @@ -78,7 +78,7 @@ void toWannier90::read_nnkp(const K_Vectors& kv) bool read_success = false; if (GlobalV::MY_RANK == 0) { - read_success = try_read_nnkp(kv); + read_success = try_read_nnkp(ucell,kv); } #ifdef __MPI @@ -87,7 +87,7 @@ void toWannier90::read_nnkp(const K_Vectors& kv) if (GlobalV::MY_RANK != 0 && read_success) { - read_success = try_read_nnkp(kv); + read_success = try_read_nnkp(ucell,kv); } #ifdef __MPI @@ -129,7 +129,7 @@ void toWannier90::cal_Mmn() { } -bool toWannier90::try_read_nnkp(const K_Vectors& kv) +bool toWannier90::try_read_nnkp(const UnitCell& ucell, const K_Vectors& kv) { std::ifstream nnkp_read(nnkpfile.c_str(), std::ios::in); @@ -144,17 +144,17 @@ bool toWannier90::try_read_nnkp(const K_Vectors& kv) nnkp_read >> real_lattice_nnkp.e11 >> real_lattice_nnkp.e12 >> real_lattice_nnkp.e13 >> real_lattice_nnkp.e21 >> real_lattice_nnkp.e22 >> real_lattice_nnkp.e23 >> real_lattice_nnkp.e31 >> real_lattice_nnkp.e32 >> real_lattice_nnkp.e33; - real_lattice_nnkp = real_lattice_nnkp / GlobalC::ucell.lat0_angstrom; - - if (std::abs(real_lattice_nnkp.e11 - GlobalC::ucell.latvec.e11) > 1.0e-4 - || std::abs(real_lattice_nnkp.e12 - GlobalC::ucell.latvec.e12) > 1.0e-4 - || std::abs(real_lattice_nnkp.e13 - GlobalC::ucell.latvec.e13) > 1.0e-4 - || std::abs(real_lattice_nnkp.e21 - GlobalC::ucell.latvec.e21) > 1.0e-4 - || std::abs(real_lattice_nnkp.e22 - GlobalC::ucell.latvec.e22) > 1.0e-4 - || std::abs(real_lattice_nnkp.e23 - GlobalC::ucell.latvec.e23) > 1.0e-4 - || std::abs(real_lattice_nnkp.e31 - GlobalC::ucell.latvec.e31) > 1.0e-4 - || std::abs(real_lattice_nnkp.e32 - GlobalC::ucell.latvec.e32) > 1.0e-4 - || std::abs(real_lattice_nnkp.e33 - GlobalC::ucell.latvec.e33) > 1.0e-4) + real_lattice_nnkp = real_lattice_nnkp / ucell.lat0_angstrom; + + if (std::abs(real_lattice_nnkp.e11 - ucell.latvec.e11) > 1.0e-4 + || std::abs(real_lattice_nnkp.e12 - ucell.latvec.e12) > 1.0e-4 + || std::abs(real_lattice_nnkp.e13 - ucell.latvec.e13) > 1.0e-4 + || std::abs(real_lattice_nnkp.e21 - ucell.latvec.e21) > 1.0e-4 + || std::abs(real_lattice_nnkp.e22 - ucell.latvec.e22) > 1.0e-4 + || std::abs(real_lattice_nnkp.e23 - ucell.latvec.e23) > 1.0e-4 + || std::abs(real_lattice_nnkp.e31 - ucell.latvec.e31) > 1.0e-4 + || std::abs(real_lattice_nnkp.e32 - ucell.latvec.e32) > 1.0e-4 + || std::abs(real_lattice_nnkp.e33 - ucell.latvec.e33) > 1.0e-4) { ModuleBase::WARNING_QUIT("toWannier90::read_nnkp", "Error real_lattice in *.nnkp file"); } @@ -170,18 +170,18 @@ bool toWannier90::try_read_nnkp(const K_Vectors& kv) nnkp_read >> recip_lattice_nnkp.e11 >> recip_lattice_nnkp.e12 >> recip_lattice_nnkp.e13 >> recip_lattice_nnkp.e21 >> recip_lattice_nnkp.e22 >> recip_lattice_nnkp.e23 >> recip_lattice_nnkp.e31 >> recip_lattice_nnkp.e32 >> recip_lattice_nnkp.e33; - const double tpiba_angstrom = ModuleBase::TWO_PI / GlobalC::ucell.lat0_angstrom; + const double tpiba_angstrom = ModuleBase::TWO_PI / ucell.lat0_angstrom; recip_lattice_nnkp = recip_lattice_nnkp / tpiba_angstrom; - if (std::abs(recip_lattice_nnkp.e11 - GlobalC::ucell.G.e11) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e12 - GlobalC::ucell.G.e12) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e13 - GlobalC::ucell.G.e13) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e21 - GlobalC::ucell.G.e21) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e22 - GlobalC::ucell.G.e22) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e23 - GlobalC::ucell.G.e23) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e31 - GlobalC::ucell.G.e31) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e32 - GlobalC::ucell.G.e32) > 1.0e-4 - || std::abs(recip_lattice_nnkp.e33 - GlobalC::ucell.G.e33) > 1.0e-4) + if (std::abs(recip_lattice_nnkp.e11 - ucell.G.e11) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e12 - ucell.G.e12) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e13 - ucell.G.e13) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e21 - ucell.G.e21) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e22 - ucell.G.e22) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e23 - ucell.G.e23) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e31 - ucell.G.e31) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e32 - ucell.G.e32) > 1.0e-4 + || std::abs(recip_lattice_nnkp.e33 - ucell.G.e33) > 1.0e-4) { ModuleBase::WARNING_QUIT("toWannier90::read_nnkp", "Error recip_lattice in *.nnkp file"); } @@ -313,7 +313,7 @@ bool toWannier90::try_read_nnkp(const K_Vectors& kv) for (int i = 0; i < num_wannier; i++) { - R_centre[i] = R_centre[i] * GlobalC::ucell.latvec; + R_centre[i] = R_centre[i] * ucell.latvec; m[i] = m[i] - 1; } diff --git a/source/module_io/to_wannier90.h b/source/module_io/to_wannier90.h index a2672aae59..a004dfada6 100644 --- a/source/module_io/to_wannier90.h +++ b/source/module_io/to_wannier90.h @@ -37,7 +37,7 @@ class toWannier90 ~toWannier90(); void calculate(); - void read_nnkp(const K_Vectors& kv); + void read_nnkp(const UnitCell& ucell, const K_Vectors& kv); void out_eig(const ModuleBase::matrix& ekb); void cal_Amn(); @@ -45,7 +45,7 @@ class toWannier90 void out_unk(); protected: - bool try_read_nnkp(const K_Vectors& kv); + bool try_read_nnkp(const UnitCell& ucell, const K_Vectors& kv); // Parameters related to k point int num_kpts; diff --git a/source/module_io/to_wannier90_lcao.cpp b/source/module_io/to_wannier90_lcao.cpp index 4d5c6052c3..72651d6991 100644 --- a/source/module_io/to_wannier90_lcao.cpp +++ b/source/module_io/to_wannier90_lcao.cpp @@ -38,14 +38,15 @@ toWannier90_LCAO::~toWannier90_LCAO() { } -void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, +void toWannier90_LCAO::calculate(const UnitCell& ucell, + const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi>& psi, const Parallel_Orbitals* pv) { this->ParaV = pv; - read_nnkp(kv); + read_nnkp(ucell,kv); if (PARAM.inp.nspin == 2) { @@ -74,11 +75,11 @@ void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, std::map>> temp_orb_index; int count = 0; - for (int it = 0; it < GlobalC::ucell.ntype; it++) + for (int it = 0; it < ucell.ntype; it++) { - for (int iL = 0; iL < GlobalC::ucell.atoms[it].nwl + 1; iL++) + for (int iL = 0; iL < ucell.atoms[it].nwl + 1; iL++) { - for (int iN = 0; iN < GlobalC::ucell.atoms[it].l_nchi[iL]; iN++) + for (int iN = 0; iN < ucell.atoms[it].l_nchi[iL]; iN++) { temp_orb_index[it][iL][iN] = count; count++; @@ -87,13 +88,13 @@ void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, } int iw = 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++) { - for (int iL = 0; iL < GlobalC::ucell.atoms[it].nwl + 1; iL++) + for (int iL = 0; iL < ucell.atoms[it].nwl + 1; iL++) { - for (int iN = 0; iN < GlobalC::ucell.atoms[it].l_nchi[iL]; iN++) + for (int iN = 0; iN < ucell.atoms[it].l_nchi[iL]; iN++) { for (int im = 0; im < (2 * iL + 1); im++) { @@ -111,10 +112,10 @@ void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, } } - initialize_orb_table(); + initialize_orb_table(ucell); produce_basis_orb(); - set_R_coor(); - count_delta_k(kv); + set_R_coor(ucell); + count_delta_k(ucell,kv); } if (out_wannier_eig) @@ -137,16 +138,16 @@ void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, return exp_idkr; }; - FR[i].set_parameters(fr_ptr[i], &GlobalC::ucell, &orb_, &GlobalC::GridD, ParaV, 140, 110); + FR[i].set_parameters(fr_ptr[i], &ucell, &orb_, &GlobalC::GridD, ParaV, 140, 110); FR[i].calculate_FR(); } - cal_Mmn(kv, psi); + cal_Mmn(ucell,kv, psi); } if (out_wannier_amn) { - cal_Amn(kv, psi); + cal_Amn(ucell,kv, psi); } if (out_wannier_unk) @@ -155,7 +156,7 @@ void toWannier90_LCAO::calculate(const ModuleBase::matrix& ekb, } } -void toWannier90_LCAO::cal_Mmn(const K_Vectors& kv, const psi::Psi>& psi) +void toWannier90_LCAO::cal_Mmn(const UnitCell& ucell, const K_Vectors& kv, const psi::Psi>& psi) { // write .mmn file std::ofstream mmn_file; @@ -180,7 +181,7 @@ void toWannier90_LCAO::cal_Mmn(const K_Vectors& kv, const psi::Psi>& psi) +void toWannier90_LCAO::cal_Amn(const UnitCell& ucell, const K_Vectors& kv, const psi::Psi>& psi) { produce_trial_in_lcao(); construct_overlap_table_project(); - cal_orbA_overlap_R(); + cal_orbA_overlap_R(ucell); // write .amn file std::ofstream Amn_file; @@ -261,7 +262,7 @@ void toWannier90_LCAO::out_unk(const psi::Psi>& psi) { } -void toWannier90_LCAO::initialize_orb_table() +void toWannier90_LCAO::initialize_orb_table(const UnitCell& ucell) { int Lmax_used = 0; int Lmax = 0; @@ -276,7 +277,7 @@ void toWannier90_LCAO::initialize_orb_table() for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb_.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb_.get_dR(); const double dk = orb_.get_dk(); @@ -302,7 +303,7 @@ void toWannier90_LCAO::initialize_orb_table() #endif } -void toWannier90_LCAO::set_R_coor() +void toWannier90_LCAO::set_R_coor(const UnitCell& ucell) { int R_minX = int(-GlobalC::GridD.getTrueCellX()); int R_minY = int(-GlobalC::GridD.getTrueCellY()); @@ -323,14 +324,14 @@ void toWannier90_LCAO::set_R_coor() for (int iz = 0; iz < R_z; iz++) { ModuleBase::Vector3 tmpR(ix + R_minX, iy + R_minY, iz + R_minZ); - R_coor_car[count] = tmpR * GlobalC::ucell.latvec; + R_coor_car[count] = tmpR * ucell.latvec; count++; } } } } -void toWannier90_LCAO::count_delta_k(const K_Vectors& kv) +void toWannier90_LCAO::count_delta_k(const UnitCell& ucell, const K_Vectors& kv) { std::set delta_k_all_tmp; for (int ik = 0; ik < cal_num_kpts; ik++) @@ -344,8 +345,8 @@ void toWannier90_LCAO::count_delta_k(const K_Vectors& kv) int cal_ikb = ikb + start_k_index; ModuleBase::Vector3 ik_car = kv.kvec_c[ik]; - ModuleBase::Vector3 ikb_car = kv.kvec_c[ikb] + G * GlobalC::ucell.G; - Abfs::Vector3_Order dk = (ikb_car - ik_car) * GlobalC::ucell.tpiba; + ModuleBase::Vector3 ikb_car = kv.kvec_c[ikb] + G * ucell.G; + Abfs::Vector3_Order dk = (ikb_car - ik_car) * ucell.tpiba; Coordinate_3D temp_dk(dk.x, dk.y, dk.z); delta_k_all_tmp.insert(temp_dk); } @@ -362,7 +363,8 @@ void toWannier90_LCAO::count_delta_k(const K_Vectors& kv) } } -void toWannier90_LCAO::unkdotkb(const K_Vectors& kv, +void toWannier90_LCAO::unkdotkb(const UnitCell& ucell, + const K_Vectors& kv, const psi::Psi>& psi_in, const int& ik, const int& ikb, @@ -379,8 +381,8 @@ void toWannier90_LCAO::unkdotkb(const K_Vectors& kv, int R_num = R_coor_car.size(); ModuleBase::Vector3 ik_car = kv.kvec_c[ik]; - ModuleBase::Vector3 ikb_car = kv.kvec_c[ikb] + G * GlobalC::ucell.G; - Abfs::Vector3_Order dk = (ikb_car - ik_car) * GlobalC::ucell.tpiba; + ModuleBase::Vector3 ikb_car = kv.kvec_c[ikb] + G * ucell.G; + Abfs::Vector3_Order dk = (ikb_car - ik_car) * ucell.tpiba; Coordinate_3D temp_dk(dk.x, dk.y, dk.z); int delta_k_index = delta_k_all_index[temp_dk]; @@ -401,7 +403,7 @@ void toWannier90_LCAO::unkdotkb(const K_Vectors& kv, auto& matrix = tmp_FR_container->get_atom_pair(iap).get_HR_values(iR); const ModuleBase::Vector3 r_index = tmp_FR_container->get_atom_pair(iap).get_R_index(iR); ModuleBase::Vector3 dR - = ModuleBase::Vector3(r_index.x, r_index.y, r_index.z) * GlobalC::ucell.latvec; + = ModuleBase::Vector3(r_index.x, r_index.y, r_index.z) * ucell.latvec; double phase = ikb_car * dR * ModuleBase::TWO_PI; std::complex kRn_phase = std::exp(ModuleBase::IMAG_UNIT * phase); for (int i = 0; i < row_size; ++i) @@ -741,7 +743,7 @@ void toWannier90_LCAO::construct_overlap_table_project() } } -void toWannier90_LCAO::cal_orbA_overlap_R() +void toWannier90_LCAO::cal_orbA_overlap_R(const UnitCell& ucell) { int row = this->ParaV->get_row_size(); int R_num = R_coor_car.size(); @@ -781,8 +783,8 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; - ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; + ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * ucell.lat0; double overlap_o = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap(orb_center, @@ -808,8 +810,8 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; - ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; + ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -838,9 +840,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -863,9 +865,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -906,8 +908,8 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; - ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; + ModuleBase::Vector3 project_orb_center = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -947,9 +949,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -972,9 +974,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -998,9 +1000,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_pz = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index] .at(1) @@ -1030,9 +1032,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -1067,9 +1069,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, @@ -1101,9 +1103,9 @@ void toWannier90_LCAO::cal_orbA_overlap_R() { ModuleBase::Vector3 R_car = R_coor_car[iR]; ModuleBase::Vector3 orb_center - = (GlobalC::ucell.atoms[it1].tau[ia1] + R_car) * GlobalC::ucell.lat0; + = (ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; ModuleBase::Vector3 project_orb_center - = R_centre[wannier_index] * GlobalC::ucell.lat0; + = R_centre[wannier_index] * ucell.lat0; double overlap_s = center2_orb11_A[iw2iorb[orb_index_row]][wannier_index].at(0).cal_overlap( orb_center, diff --git a/source/module_io/to_wannier90_lcao.h b/source/module_io/to_wannier90_lcao.h index b33dfae9d7..34d0e64dbb 100644 --- a/source/module_io/to_wannier90_lcao.h +++ b/source/module_io/to_wannier90_lcao.h @@ -79,12 +79,14 @@ class toWannier90_LCAO : public toWannier90 ); ~toWannier90_LCAO(); - void calculate(const ModuleBase::matrix& ekb, + void calculate(const UnitCell& ucell, + const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi>& psi, const Parallel_Orbitals* pv); - void calculate(const ModuleBase::matrix& ekb, + void calculate(const UnitCell& ucell, + const ModuleBase::matrix& ekb, const K_Vectors& kv, const psi::Psi& psi, const Parallel_Orbitals* pv) @@ -92,8 +94,8 @@ class toWannier90_LCAO : public toWannier90 throw std::logic_error("The wave function of toWannier90_LCAO_IN_PW is generally a std::complex type."); } - void cal_Amn(const K_Vectors& kv, const psi::Psi>& psi); - void cal_Mmn(const K_Vectors& kv, const psi::Psi>& psi); + void cal_Amn(const UnitCell& ucell, const K_Vectors& kv, const psi::Psi>& psi); + void cal_Mmn(const UnitCell& ucell, const K_Vectors& kv, const psi::Psi>& psi); void out_unk(const psi::Psi>& psi); protected: @@ -127,15 +129,16 @@ class toWannier90_LCAO : public toWannier90 const Parallel_Orbitals* ParaV; - void initialize_orb_table(); + void initialize_orb_table(const UnitCell& ucell); void produce_basis_orb(); - void set_R_coor(); - void count_delta_k(const K_Vectors& kv); + void set_R_coor(const UnitCell& ucell); + void count_delta_k(const UnitCell& ucell, const K_Vectors& kv); std::vector delta_k_all; std::map delta_k_all_index; - void unkdotkb(const K_Vectors& kv, + void unkdotkb(const UnitCell& ucell, + const K_Vectors& kv, const psi::Psi>& psi_in, const int& ik, const int& ikb, @@ -144,7 +147,7 @@ class toWannier90_LCAO : public toWannier90 void produce_trial_in_lcao(); void construct_overlap_table_project(); - void cal_orbA_overlap_R(); + void cal_orbA_overlap_R(const UnitCell& ucell); void unkdotA(const K_Vectors& kv, const psi::Psi>& psi_in, diff --git a/source/module_io/to_wannier90_lcao_in_pw.cpp b/source/module_io/to_wannier90_lcao_in_pw.cpp index 0419009d17..bf6bac22f6 100644 --- a/source/module_io/to_wannier90_lcao_in_pw.cpp +++ b/source/module_io/to_wannier90_lcao_in_pw.cpp @@ -30,6 +30,7 @@ toWannier90_LCAO_IN_PW::~toWannier90_LCAO_IN_PW() } void toWannier90_LCAO_IN_PW::calculate( + UnitCell& ucell, const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, @@ -45,13 +46,13 @@ void toWannier90_LCAO_IN_PW::calculate( ModulePW::PW_Basis_K* wfcpw_ptr = const_cast(wfcpw); this->psi_init_ = new psi_initializer_nao, base_device::DEVICE_CPU>(); #ifdef __MPI - this->psi_init_->initialize(sf_ptr, wfcpw_ptr, &(GlobalC::ucell), &(GlobalC::Pkpoints), 1, nullptr, GlobalV::MY_RANK); + this->psi_init_->initialize(sf_ptr, wfcpw_ptr, &ucell, &(GlobalC::Pkpoints), 1, nullptr, GlobalV::MY_RANK); #else - this->psi_init_->initialize(sf_ptr, wfcpw_ptr, &(GlobalC::ucell), 1, nullptr); + this->psi_init_->initialize(sf_ptr, wfcpw_ptr, &(ucell), 1, nullptr); #endif this->psi_init_->tabulate(); this->psi_init_->allocate(true); - read_nnkp(kv); + read_nnkp(ucell,kv); if (PARAM.inp.nspin == 2) { @@ -69,7 +70,7 @@ void toWannier90_LCAO_IN_PW::calculate( } } - psi::Psi> *unk_inLcao = get_unk_from_lcao(*psi, wfcpw, sf, kv); + psi::Psi> *unk_inLcao = get_unk_from_lcao(ucell,*psi, wfcpw, sf, kv); if (out_wannier_eig) { @@ -106,6 +107,7 @@ void toWannier90_LCAO_IN_PW::calculate( } psi::Psi>* toWannier90_LCAO_IN_PW::get_unk_from_lcao( + const UnitCell& ucell, const psi::Psi>& psi_in, const ModulePW::PW_Basis_K* wfcpw, const Structure_Factor& sf, @@ -118,7 +120,7 @@ psi::Psi>* toWannier90_LCAO_IN_PW::get_unk_from_lcao( unk_inLcao->zero_out(); // Orbital projection to plane wave - ModuleBase::realArray table_local(GlobalC::ucell.ntype, GlobalC::ucell.nmax_total, PARAM.globalv.nqx); + ModuleBase::realArray table_local(ucell.ntype, ucell.nmax_total, PARAM.globalv.nqx); for (int ik = 0; ik < num_kpts; ik++) { diff --git a/source/module_io/to_wannier90_lcao_in_pw.h b/source/module_io/to_wannier90_lcao_in_pw.h index 5362fc3792..26fe2abb92 100644 --- a/source/module_io/to_wannier90_lcao_in_pw.h +++ b/source/module_io/to_wannier90_lcao_in_pw.h @@ -46,7 +46,8 @@ class toWannier90_LCAO_IN_PW : public toWannier90_PW const std::string& wannier_spin); ~toWannier90_LCAO_IN_PW(); - void calculate(const ModuleBase::matrix& ekb, + void calculate(UnitCell& ucell, + const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, const Structure_Factor& sf, @@ -54,7 +55,8 @@ class toWannier90_LCAO_IN_PW : public toWannier90_PW const psi::Psi>* psi, const Parallel_Orbitals* pv); - void calculate(const ModuleBase::matrix& ekb, + void calculate(UnitCell& ucell, + const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, const Structure_Factor& sf, @@ -76,7 +78,8 @@ class toWannier90_LCAO_IN_PW : public toWannier90_PW /// @param sf [in] computational methods instance, structure factor calculator /// @param kv [in] data carrier, storing kpoints information /// @return psi::Psi>* - psi::Psi>* get_unk_from_lcao(const psi::Psi>& psi_in, + psi::Psi>* get_unk_from_lcao(const UnitCell& ucell, + const psi::Psi>& psi_in, const ModulePW::PW_Basis_K* wfcpw, const Structure_Factor& sf, const K_Vectors& kv); diff --git a/source/module_io/to_wannier90_pw.cpp b/source/module_io/to_wannier90_pw.cpp index 9afb105d9d..e1ab2477f2 100644 --- a/source/module_io/to_wannier90_pw.cpp +++ b/source/module_io/to_wannier90_pw.cpp @@ -28,6 +28,7 @@ toWannier90_PW::~toWannier90_PW() } void toWannier90_PW::calculate( + const UnitCell& ucell, const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, @@ -35,7 +36,7 @@ void toWannier90_PW::calculate( const psi::Psi>* psi ) { - read_nnkp(kv); + read_nnkp(ucell,kv); if (PARAM.inp.nspin == 2) { @@ -73,6 +74,12 @@ void toWannier90_PW::calculate( out_unk(*psi, wfcpw, bigpw); } +} +void toWannier90_PW::set_tpiba_omega(const double& tpiba, const double& omega) +{ + this->tpiba = &tpiba; + this->omega = ω + } void toWannier90_PW::cal_Mmn( @@ -946,9 +953,10 @@ void toWannier90_PW::get_trial_orbitals_lm_k( std::complex *orbital_in_G_single ) { + const double tpiba = *this->tpiba; for (int ig = 0; ig < npw; ig++) { - orbital_in_G_single[ig] = ModuleBase::PolyInt::Polynomial_Interpolation(radial_in_q_single, PARAM.globalv.nqx, PARAM.globalv.dq, gk[ig].norm() * GlobalC::ucell.tpiba); + orbital_in_G_single[ig] = ModuleBase::PolyInt::Polynomial_Interpolation(radial_in_q_single, PARAM.globalv.nqx, PARAM.globalv.dq, gk[ig].norm() * tpiba); } std::complex lphase = pow(ModuleBase::NEG_IMAG_UNIT, orbital_L); @@ -970,7 +978,7 @@ void toWannier90_PW::integral( double *table ) { - const double pref = ModuleBase::FOUR_PI / sqrt(GlobalC::ucell.omega); + const double pref = ModuleBase::FOUR_PI / sqrt(*this->omega); double *inner_part = new double[meshr]; for (int ir = 0; ir < meshr; ir++) diff --git a/source/module_io/to_wannier90_pw.h b/source/module_io/to_wannier90_pw.h index d6ab7ccc12..f7f5768ff2 100644 --- a/source/module_io/to_wannier90_pw.h +++ b/source/module_io/to_wannier90_pw.h @@ -34,6 +34,7 @@ class toWannier90_PW : public toWannier90 ~toWannier90_PW(); void calculate( + const UnitCell& ucell, const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, @@ -42,6 +43,7 @@ class toWannier90_PW : public toWannier90 ); void calculate( + const UnitCell& ucell, const ModuleBase::matrix& ekb, const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw, @@ -59,12 +61,14 @@ class toWannier90_PW : public toWannier90 const ModulePW::PW_Basis_K* wfcpw, const ModulePW::PW_Basis_Big* bigpw ); - + void set_tpiba_omega(const double& tpiba, const double& omega); protected: // Radial section of trial orbitals const int mesh_r = 333; const double dx = 0.025; const double x_min = -6.0; + double const *tpiba; + double const *omega; void unkdotkb( const psi::Psi>& psi_pw, diff --git a/source/module_io/unk_overlap_lcao.cpp b/source/module_io/unk_overlap_lcao.cpp index b8ca83ebf2..b784f54a1c 100644 --- a/source/module_io/unk_overlap_lcao.cpp +++ b/source/module_io/unk_overlap_lcao.cpp @@ -25,7 +25,10 @@ unkOverlap_lcao::~unkOverlap_lcao() // GlobalV::ofs_running << "this is ~unkOverlap_lcao()" << std::endl; } -void unkOverlap_lcao::init(const Grid_Technique& gt, const int nkstot, const LCAO_Orbitals& orb) +void unkOverlap_lcao::init(const UnitCell& ucell, + const Grid_Technique& gt, + const int nkstot, + const LCAO_Orbitals& orb) { // std::cout << "unkOverlap_lcao::init start" << std::endl; @@ -40,7 +43,7 @@ void unkOverlap_lcao::init(const Grid_Technique& gt, const int nkstot, const LCA for (int it = 0; it < ntype; it++) { lmax_orb = std::max(lmax_orb, orb.Phi[it].getLmax()); - lmax_beta = std::max(lmax_beta, GlobalC::ucell.infoNL.Beta[it].getLmax()); + lmax_beta = std::max(lmax_beta, ucell.infoNL.Beta[it].getLmax()); } const double dr = orb.get_dR(); const double dk = orb.get_dk(); @@ -140,9 +143,9 @@ void unkOverlap_lcao::init(const Grid_Technique& gt, const int nkstot, const LCA } } - for (int TA = 0; TA < GlobalC::ucell.ntype; TA++) + for (int TA = 0; TA < ucell.ntype; TA++) { - for (int TB = 0; TB < GlobalC::ucell.ntype; TB++) + for (int TB = 0; TB < ucell.ntype; TB++) { for (int LA = 0; LA <= orb.Phi[TA].getLmax(); LA++) { @@ -165,9 +168,9 @@ void unkOverlap_lcao::init(const Grid_Technique& gt, const int nkstot, const LCA } } - for (int TA = 0; TA < GlobalC::ucell.ntype; TA++) + for (int TA = 0; TA < ucell.ntype; TA++) { - for (int TB = 0; TB < GlobalC::ucell.ntype; TB++) + for (int TB = 0; TB < ucell.ntype; TB++) { for (int LA = 0; LA <= orb.Phi[TA].getLmax(); LA++) { @@ -216,18 +219,18 @@ void unkOverlap_lcao::init(const Grid_Technique& gt, const int nkstot, const LCA return; } -int unkOverlap_lcao::iw2it(int iw) +int unkOverlap_lcao::iw2it(const UnitCell& ucell, int iw) { int ic, type; ic = 0; type = 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++) { - for (int L = 0; L < GlobalC::ucell.atoms[it].nwl + 1; L++) + for (int L = 0; L < ucell.atoms[it].nwl + 1; L++) { - 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++) { for (int i = 0; i < (2 * L + 1); i++) { @@ -244,18 +247,18 @@ int unkOverlap_lcao::iw2it(int iw) return type; } -int unkOverlap_lcao::iw2ia(int iw) +int unkOverlap_lcao::iw2ia(const UnitCell& ucell,int iw) { int ic, na; ic = 0; na = 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++) { - for (int L = 0; L < GlobalC::ucell.atoms[it].nwl + 1; L++) + for (int L = 0; L < ucell.atoms[it].nwl + 1; L++) { - 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++) { for (int i = 0; i < (2 * L + 1); i++) { @@ -272,18 +275,18 @@ int unkOverlap_lcao::iw2ia(int iw) return na; } -int unkOverlap_lcao::iw2iL(int iw) +int unkOverlap_lcao::iw2iL(const UnitCell& ucell, int iw) { int ic, iL; ic = 0; iL = 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++) { - for (int L = 0; L < GlobalC::ucell.atoms[it].nwl + 1; L++) + for (int L = 0; L < ucell.atoms[it].nwl + 1; L++) { - 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++) { for (int i = 0; i < (2 * L + 1); i++) { @@ -300,18 +303,18 @@ int unkOverlap_lcao::iw2iL(int iw) return iL; } -int unkOverlap_lcao::iw2iN(int iw) +int unkOverlap_lcao::iw2iN(const UnitCell& ucell,int iw) { int ic, iN; ic = 0; iN = 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++) { - for (int L = 0; L < GlobalC::ucell.atoms[it].nwl + 1; L++) + for (int L = 0; L < ucell.atoms[it].nwl + 1; L++) { - 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++) { for (int i = 0; i < (2 * L + 1); i++) { @@ -328,18 +331,18 @@ int unkOverlap_lcao::iw2iN(int iw) return iN; } -int unkOverlap_lcao::iw2im(int iw) +int unkOverlap_lcao::iw2im(const UnitCell& ucell, int iw) { int ic, im; ic = 0; im = 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++) { - for (int L = 0; L < GlobalC::ucell.atoms[it].nwl + 1; L++) + for (int L = 0; L < ucell.atoms[it].nwl + 1; L++) { - 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++) { for (int i = 0; i < (2 * L + 1); i++) { @@ -357,7 +360,7 @@ int unkOverlap_lcao::iw2im(int iw) } // search for the nearest neighbor atoms -void unkOverlap_lcao::cal_R_number() +void unkOverlap_lcao::cal_R_number(const UnitCell& ucell) { // The number of overlaps between atomic orbitals 1 and atomic orbitals 2, // or the number of R, is empty when there is no overlap @@ -368,39 +371,39 @@ void unkOverlap_lcao::cal_R_number() } ModuleBase::Vector3 tau1, tau2, dtau; - 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(GlobalC::ucell, tau1, T1, I1); + GlobalC::GridD.Find_atom(ucell, tau1, T1, I1); for (int ad = 0; ad < GlobalC::GridD.getAdjacentNum() + 1; ++ad) { const int T2 = GlobalC::GridD.getType(ad); const int I2 = GlobalC::GridD.getNatom(ad); - Atom* atom2 = &GlobalC::ucell.atoms[T2]; + Atom* atom2 = &ucell.atoms[T2]; const double R_direct_x = (double)GlobalC::GridD.getBox(ad).x; const double R_direct_y = (double)GlobalC::GridD.getBox(ad).y; const double R_direct_z = (double)GlobalC::GridD.getBox(ad).z; tau2 = GlobalC::GridD.getAdjacentTau(ad); dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::ucell.lat0; + double distance = dtau.norm() * ucell.lat0; double rcut = rcut_orb_[T1] + rcut_orb_[T2]; if (distance < rcut - 1.0e-15) { - // translate: the unit of R_car is GlobalC::ucell.lat0 - ModuleBase::Vector3 R_car = R_direct_x * GlobalC::ucell.a1 + R_direct_y * GlobalC::ucell.a2 - + R_direct_z * GlobalC::ucell.a3; + // translate: the unit of R_car is ucell.lat0 + ModuleBase::Vector3 R_car = R_direct_x * ucell.a1 + R_direct_y * ucell.a2 + + R_direct_z * ucell.a3; for (int iw1 = 0; iw1 < atom1->nw; iw1++) { - int orb_index_in_NLOCAL_1 = GlobalC::ucell.itiaiw2iwt(T1, I1, iw1); + int orb_index_in_NLOCAL_1 = ucell.itiaiw2iwt(T1, I1, iw1); for (int iw2 = 0; iw2 < atom2->nw; iw2++) { - int orb_index_in_NLOCAL_2 = GlobalC::ucell.itiaiw2iwt(T2, I2, iw2); + int orb_index_in_NLOCAL_2 = ucell.itiaiw2iwt(T2, I2, iw2); orb1_orb2_R[orb_index_in_NLOCAL_1][orb_index_in_NLOCAL_2].push_back(R_car); } // end iw2 @@ -416,7 +419,7 @@ void unkOverlap_lcao::cal_R_number() return; } -void unkOverlap_lcao::cal_orb_overlap() +void unkOverlap_lcao::cal_orb_overlap(const UnitCell& ucell) { // std::cout << "the cal_orb_overlap is start" << std::endl; psi_psi.resize(PARAM.globalv.nlocal); @@ -439,23 +442,23 @@ void unkOverlap_lcao::cal_orb_overlap() if (orb1_orb2_R[iw1][iw2].empty()) continue; - int atomType1 = iw2it(iw1); - int ia1 = iw2ia(iw1); - int N1 = iw2iN(iw1); - int L1 = iw2iL(iw1); - int m1 = iw2im(iw1); - int atomType2 = iw2it(iw2); - int ia2 = iw2ia(iw2); - int N2 = iw2iN(iw2); - int L2 = iw2iL(iw2); - int m2 = iw2im(iw2); + int atomType1 = iw2it(ucell,iw1); + int ia1 = iw2ia(ucell,iw1); + int N1 = iw2iN(ucell,iw1); + int L1 = iw2iL(ucell,iw1); + int m1 = iw2im(ucell,iw1); + int atomType2 = iw2it(ucell,iw2); + int ia2 = iw2ia(ucell,iw2); + int N2 = iw2iN(ucell,iw2); + int L2 = iw2iL(ucell,iw2); + int m2 = iw2im(ucell,iw2); for (int iR = 0; iR < orb1_orb2_R[iw1][iw2].size(); iR++) { ModuleBase::Vector3 r_distance - = (GlobalC::ucell.atoms[atomType2].tau[ia2] - GlobalC::ucell.atoms[atomType1].tau[ia1] + = (ucell.atoms[atomType2].tau[ia2] - ucell.atoms[atomType1].tau[ia1] + orb1_orb2_R[iw1][iw2][iR]) - * GlobalC::ucell.lat0; + * ucell.lat0; psi_psi[iw1][iw2].push_back( center2_orb11[atomType1][atomType2][L1][N1][L2].at(N2).cal_overlap(origin_point, r_distance, @@ -491,7 +494,8 @@ void unkOverlap_lcao::cal_orb_overlap() return; } -void unkOverlap_lcao::prepare_midmatrix_pblas(const int ik_L, +void unkOverlap_lcao::prepare_midmatrix_pblas(const UnitCell& ucell, + const int ik_L, const int ik_R, const ModuleBase::Vector3 dk, std::complex*& midmatrix, @@ -511,13 +515,13 @@ void unkOverlap_lcao::prepare_midmatrix_pblas(const int ik_L, if (ir >= 0 && ic >= 0) { int index = ic * pv.nrow + ir; - ModuleBase::Vector3 tau1 = GlobalC::ucell.atoms[iw2it(iw_row)].tau[iw2ia(iw_row)]; + ModuleBase::Vector3 tau1 = ucell.atoms[iw2it(ucell,iw_row)].tau[iw2ia(ucell,iw_row)]; for (int iR = 0; iR < orb1_orb2_R[iw_row][iw_col].size(); iR++) { double kRn = (kv.kvec_c[ik_R] * orb1_orb2_R[iw_row][iw_col][iR] - dk * tau1) * ModuleBase::TWO_PI; std::complex kRn_phase(cos(kRn), sin(kRn)); std::complex orb_overlap(psi_psi[iw_row][iw_col][iR], - (-dk * GlobalC::ucell.tpiba * psi_r_psi[iw_row][iw_col][iR])); + (-dk * ucell.tpiba * psi_r_psi[iw_row][iw_col][iR])); midmatrix[index] = midmatrix[index] + kRn_phase * orb_overlap; } } @@ -525,7 +529,8 @@ void unkOverlap_lcao::prepare_midmatrix_pblas(const int ik_L, } } -std::complex unkOverlap_lcao::det_berryphase(const int ik_L, +std::complex unkOverlap_lcao::det_berryphase(const UnitCell& ucell, + const int ik_L, const int ik_R, const ModuleBase::Vector3 dk, const int occ_bands, @@ -541,7 +546,7 @@ std::complex unkOverlap_lcao::det_berryphase(const int ik_L, ModuleBase::GlobalFunc::ZEROS(C_matrix, para_orb.nloc); ModuleBase::GlobalFunc::ZEROS(out_matrix, para_orb.nloc); - this->prepare_midmatrix_pblas(ik_L, ik_R, dk, midmatrix, para_orb, kv); + this->prepare_midmatrix_pblas(ucell,ik_L, ik_R, dk, midmatrix, para_orb, kv); char transa = 'C'; char transb = 'N'; diff --git a/source/module_io/unk_overlap_lcao.h b/source/module_io/unk_overlap_lcao.h index 8f5aad7cf8..0ba664271e 100644 --- a/source/module_io/unk_overlap_lcao.h +++ b/source/module_io/unk_overlap_lcao.h @@ -48,21 +48,23 @@ class unkOverlap_lcao unkOverlap_lcao(); ~unkOverlap_lcao(); - void init(const Grid_Technique& gt, const int nkstot, const LCAO_Orbitals& orb); - int iw2it(int iw); - int iw2ia(int iw); - int iw2iL(int iw); - int iw2iN(int iw); - int iw2im(int iw); - void cal_R_number(); - void cal_orb_overlap(); - void prepare_midmatrix_pblas(const int ik_L, + void init(const UnitCell& ucell, const Grid_Technique& gt, const int nkstot, const LCAO_Orbitals& orb); + int iw2it(const UnitCell& ucell, int iw); + int iw2ia(const UnitCell& ucell, int iw); + int iw2iL(const UnitCell& ucell, int iw); + int iw2iN(const UnitCell& ucell, int iw); + int iw2im(const UnitCell& ucell, int iw); + void cal_R_number(const UnitCell& ucell); + void cal_orb_overlap(const UnitCell& ucell); + void prepare_midmatrix_pblas(const UnitCell& ucell, + const int ik_L, const int ik_R, const ModuleBase::Vector3 dk, std::complex*& midmatrix, const Parallel_Orbitals& pv, const K_Vectors& kv); - std::complex det_berryphase(const int ik_L, + std::complex det_berryphase(const UnitCell& ucell, + const int ik_L, const int ik_R, const ModuleBase::Vector3 dk, const int occ_bands, diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index 57794e45ed..e1a7c0fa4d 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -5,7 +5,8 @@ #include "module_io/dipole_io.h" // fuxiang add 2017-03-15 -void ModuleIO::write_dipole(const double* rho_save, +void ModuleIO::write_dipole(const UnitCell& ucell, + const double* rho_save, const ModulePW::PW_Basis* rhopw, const int& is, const int& istep, @@ -32,14 +33,14 @@ void ModuleIO::write_dipole(const double* rho_save, double bmod[3]; for (int i = 0; i < 3; i++) { - bmod[i] = prepare(GlobalC::ucell, i); + bmod[i] = prepare(ucell, i); } #ifndef __MPI double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; - double lat_factor_x = GlobalC::ucell.lat0 * 0.529177 / rhopw->nx; - double lat_factor_y = GlobalC::ucell.lat0 * 0.529177 / rhopw->ny; - double lat_factor_z = GlobalC::ucell.lat0 * 0.529177 / rhopw->nz; + double lat_factor_x = ucell.lat0 * 0.529177 / rhopw->nx; + double lat_factor_y = ucell.lat0 * 0.529177 / rhopw->ny; + double lat_factor_z = ucell.lat0 * 0.529177 / rhopw->nz; for (int k = 0; k < rhopw->nz; k++) { @@ -57,9 +58,9 @@ void ModuleIO::write_dipole(const double* rho_save, } } - dipole_elec_x *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); - dipole_elec_y *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); - dipole_elec_z *= GlobalC::ucell.omega / static_cast(rhopw->nxyz); + dipole_elec_x *= ucell.omega / static_cast(rhopw->nxyz); + dipole_elec_y *= ucell.omega / static_cast(rhopw->nxyz); + dipole_elec_z *= ucell.omega / static_cast(rhopw->nxyz); Parallel_Reduce::reduce_pool(dipole_elec_x); Parallel_Reduce::reduce_pool(dipole_elec_y); Parallel_Reduce::reduce_pool(dipole_elec_z); @@ -88,7 +89,7 @@ void ModuleIO::write_dipole(const double* rho_save, Parallel_Reduce::reduce_pool(dipole_elec[2]); for (int i = 0; i < 3; ++i) { - dipole_elec[i] *= GlobalC::ucell.lat0 / bmod[i] * GlobalC::ucell.omega / rhopw->nxyz; + dipole_elec[i] *= ucell.lat0 / bmod[i] * ucell.omega / rhopw->nxyz; } std::cout << std::setprecision(15) << "dipole_elec_x: " << dipole_elec[0] << std::endl; @@ -103,16 +104,16 @@ void ModuleIO::write_dipole(const double* rho_save, for (int i = 0; i < 3; ++i) { - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < ucell.ntype; ++it) { double sum = 0; - for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ++ia) + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) { - sum += GlobalC::ucell.atoms[it].taud[ia][i]; + sum += ucell.atoms[it].taud[ia][i]; } - dipole_ion[i] += sum * GlobalC::ucell.atoms[it].ncpp.zv; + dipole_ion[i] += sum * ucell.atoms[it].ncpp.zv; } - dipole_ion[i] *= GlobalC::ucell.lat0 / bmod[i]; //* ModuleBase::FOUR_PI / GlobalC::ucell.omega; + dipole_ion[i] *= ucell.lat0 / bmod[i]; //* ModuleBase::FOUR_PI / ucell.omega; } std::cout << std::setprecision(8) << "dipole_ion_x: " << dipole_ion[0] << std::endl; diff --git a/source/module_io/write_dmr.cpp b/source/module_io/write_dmr.cpp index 28b5c00a25..b7ce1124eb 100644 --- a/source/module_io/write_dmr.cpp +++ b/source/module_io/write_dmr.cpp @@ -55,6 +55,8 @@ void write_dmr(const std::vector*> dmr, const bool out_csr, const bool out_npz, const bool append, + const int* iat2iwt, + const int* nat, const int istep) { if (!out_csr && !out_npz) @@ -73,7 +75,7 @@ void write_dmr(const std::vector*> dmr, Parallel_Orbitals serialV; serialV.init(nbasis, nbasis, nbasis, paraV.comm()); serialV.set_serial(nbasis, nbasis); - serialV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), GlobalC::ucell.nat, nbasis); + serialV.set_atomic_trace(iat2iwt, *nat, nbasis); hamilt::HContainer dm_serial(&serialV); hamilt::gatherParallels(*dmr[ispin], &dm_serial, 0); #else diff --git a/source/module_io/write_dmr.h b/source/module_io/write_dmr.h index d3cc1b4573..24077bffca 100644 --- a/source/module_io/write_dmr.h +++ b/source/module_io/write_dmr.h @@ -45,6 +45,8 @@ void write_dmr(const std::vector*> dmr, const bool out_csr, const bool out_npz, const bool append, + const int* iat2iwt, + const int* nat, const int istep); } // namespace ModuleIO diff --git a/source/module_io/write_dos_lcao.cpp b/source/module_io/write_dos_lcao.cpp index 89dda2e591..e475c77459 100644 --- a/source/module_io/write_dos_lcao.cpp +++ b/source/module_io/write_dos_lcao.cpp @@ -27,7 +27,8 @@ #include "module_base/scalapack_connector.h" template <> -void ModuleIO::write_dos_lcao(const psi::Psi* psi, +void ModuleIO::write_dos_lcao(const UnitCell& ucell, + const psi::Psi* psi, const Parallel_Orbitals& pv, const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg, @@ -254,24 +255,24 @@ void ModuleIO::write_dos_lcao(const psi::Psi* psi, out << std::setw(20) << en << std::endl; } out << "" << std::endl; - for (int i = 0; i < GlobalC::ucell.nat; i++) + for (int i = 0; i < ucell.nat; i++) { - int a = GlobalC::ucell.iat2ia[i]; - int t = GlobalC::ucell.iat2it[i]; - Atom* atom1 = &GlobalC::ucell.atoms[t]; - const int s0 = GlobalC::ucell.itiaiw2iwt(t, a, 0); + int a = ucell.iat2ia[i]; + int t = ucell.iat2it[i]; + Atom* atom1 = &ucell.atoms[t]; + const int s0 = ucell.itiaiw2iwt(t, a, 0); for (int j = 0; j < atom1->nw; ++j) { const int L1 = atom1->iw2l[j]; const int N1 = atom1->iw2n[j]; const int m1 = atom1->iw2m[j]; - const int w = GlobalC::ucell.itiaiw2iwt(t, a, j); + const int w = ucell.itiaiw2iwt(t, a, j); // out << "" <* psi, out << "" << std::endl; out.close(); } - ModuleIO::write_orb_info(&(GlobalC::ucell)); + ModuleIO::write_orb_info(&(ucell)); } delete[] pdos; @@ -341,7 +342,8 @@ void ModuleIO::write_dos_lcao(const psi::Psi* psi, } template <> -void ModuleIO::write_dos_lcao(const psi::Psi>* psi, +void ModuleIO::write_dos_lcao(const UnitCell& ucell, + const psi::Psi>* psi, const Parallel_Orbitals& pv, const ModuleBase::matrix& ekb, const ModuleBase::matrix& wg, @@ -604,24 +606,24 @@ void ModuleIO::write_dos_lcao(const psi::Psi>* psi, out << std::setw(20) << en << std::endl; } out << "" << std::endl; - for (int i = 0; i < GlobalC::ucell.nat; i++) + for (int i = 0; i < ucell.nat; i++) { - int a = GlobalC::ucell.iat2ia[i]; - int t = GlobalC::ucell.iat2it[i]; - Atom* atom1 = &GlobalC::ucell.atoms[t]; - const int s0 = GlobalC::ucell.itiaiw2iwt(t, a, 0); + int a = ucell.iat2ia[i]; + int t = ucell.iat2it[i]; + Atom* atom1 = &ucell.atoms[t]; + const int s0 = ucell.itiaiw2iwt(t, a, 0); for (int j = 0; j < atom1->nw; ++j) { const int L1 = atom1->iw2l[j]; const int N1 = atom1->iw2n[j]; const int m1 = atom1->iw2m[j]; - const int w = GlobalC::ucell.itiaiw2iwt(t, a, j); + const int w = ucell.itiaiw2iwt(t, a, j); // out << "" <>* psi, out << "" << std::endl; out.close(); } - ModuleIO::write_orb_info(&(GlobalC::ucell)); + ModuleIO::write_orb_info(&(ucell)); } delete[] pdos; } diff --git a/source/module_io/write_dos_lcao.h b/source/module_io/write_dos_lcao.h index ebfd057a34..ec4e398b2a 100644 --- a/source/module_io/write_dos_lcao.h +++ b/source/module_io/write_dos_lcao.h @@ -12,6 +12,7 @@ namespace ModuleIO /// @brief calculate density of states(DOS) and partial density of states(PDOS) and mulliken charge for LCAO base template void write_dos_lcao( + const UnitCell& ucell, const psi::Psi* psi, const Parallel_Orbitals &pv, const ModuleBase::matrix& ekb, diff --git a/source/module_io/write_elecstat_pot.cpp b/source/module_io/write_elecstat_pot.cpp index eaf7762285..f4561bb606 100644 --- a/source/module_io/write_elecstat_pot.cpp +++ b/source/module_io/write_elecstat_pot.cpp @@ -100,7 +100,7 @@ void write_elecstat_pot( istep, fn, ef_tmp, - &(GlobalC::ucell), + ucell, precision, out_fermi); diff --git a/source/module_io/write_vxc_lip.hpp b/source/module_io/write_vxc_lip.hpp index 94890902d0..f3c30bd5e2 100644 --- a/source/module_io/write_vxc_lip.hpp +++ b/source/module_io/write_vxc_lip.hpp @@ -197,7 +197,7 @@ namespace ModuleIO // for (int ir = 0;ir < potxc->get_veff_smooth().nc;++ir) // exc_by_rho += potxc->get_veff_smooth()(0, ir) * chg.rho[0][ir]; // Parallel_Reduce::reduce_all(exc_by_rho); - // exc_by_rho *= ((FPTYPE)GlobalC::ucell.omega * (FPTYPE)GlobalV::NPROC / (FPTYPE)potxc->get_veff_smooth().nc); + // exc_by_rho *= ((FPTYPE)ucell.omega * (FPTYPE)GlobalV::NPROC / (FPTYPE)potxc->get_veff_smooth().nc); // std::cout << "xc all-bands energy by rho =" << exc_by_rho << std::endl; //===== test total xc energy ======= //===== test total exx energy ======= diff --git a/source/module_io/write_wfc_r.cpp b/source/module_io/write_wfc_r.cpp index b900f8f0be..4b171f26f6 100644 --- a/source/module_io/write_wfc_r.cpp +++ b/source/module_io/write_wfc_r.cpp @@ -23,7 +23,8 @@ namespace ModuleIO // write ||wfc_r|| for all k-points and all bands // Input: wfc_g(ik, ib, ig) // loop order is for(z){for(y){for(x)}} -void write_psi_r_1(const psi::Psi>& wfc_g, +void write_psi_r_1(const UnitCell& ucell, + const psi::Psi>& wfc_g, const ModulePW::PW_Basis_K* wfcpw, const std::string& folder_name, const bool& square, @@ -66,18 +67,16 @@ void write_psi_r_1(const psi::Psi>& wfc_g, wfc_imag[ir] = wfc_r[ir].imag(); } } - const std::string file_name = outdir + "wfc_realspace_" + ModuleBase::GlobalFunc::TO_STRING(ik_out) + "_" - + ModuleBase::GlobalFunc::TO_STRING(ib); + const std::string file_name_base = outdir + "wfc_realspace_" + + ModuleBase::GlobalFunc::TO_STRING(ik_out) + "_" + + ModuleBase::GlobalFunc::TO_STRING(ib); + const std::string file_name = square ? file_name_base : file_name_base + "_imag"; #ifdef __MPI // Use write_chg_r_1 to output the real and imaginary parts of the wave function to file mpi_requests.push_back({}); - write_chg_r_1(wfcpw, wfc_real, file_name, mpi_requests.back()); - if (!square) - { - write_chg_r_1(wfcpw, wfc_imag, file_name + "_imag", mpi_requests.back()); - } + write_chg_r_1(ucell,wfcpw, wfc_real, file_name, mpi_requests.back()); #else - write_chg_r_1(wfcpw, wfc_real, file_name); + write_chg_r_1(ucell,wfcpw, wfc_real, file_name); // if (!square) // write_chg_r_1(wfc_imag, file_name + "_imag", mpi_requests.back()); #endif @@ -115,14 +114,14 @@ std::vector> cal_wfc_r(const ModulePW::PW_Basis_K* wfcpw, } // Input: chg_r[ir] -#ifdef __MPI -void write_chg_r_1(const ModulePW::PW_Basis_K* wfcpw, +void write_chg_r_1(const UnitCell& ucell, + const ModulePW::PW_Basis_K* wfcpw, const std::vector& chg_r, - const std::string& file_name, - MPI_Request& mpi_request) -#else -void write_chg_r_1(const ModulePW::PW_Basis_K* wfcpw, const std::vector& chg_r, const std::string& file_name) -#endif + const std::string& file_name + #ifdef __MPI + ,MPI_Request& mpi_request + #endif + ) { ModuleBase::timer::tick("ModuleIO", "write_chg_r_1"); std::ofstream ofs; @@ -135,32 +134,32 @@ void write_chg_r_1(const ModulePW::PW_Basis_K* wfcpw, const std::vector& ofs.open(file_name); ofs << "calculated by ABACUS" << std::endl; - ofs << GlobalC::ucell.lat0_angstrom << std::endl; - ofs << GlobalC::ucell.latvec.e11 << " " << GlobalC::ucell.latvec.e12 << " " << GlobalC::ucell.latvec.e13 + ofs << ucell.lat0_angstrom << std::endl; + ofs << ucell.latvec.e11 << " " << ucell.latvec.e12 << " " << ucell.latvec.e13 << std::endl - << GlobalC::ucell.latvec.e21 << " " << GlobalC::ucell.latvec.e22 << " " << GlobalC::ucell.latvec.e23 + << ucell.latvec.e21 << " " << ucell.latvec.e22 << " " << ucell.latvec.e23 << std::endl - << GlobalC::ucell.latvec.e31 << " " << GlobalC::ucell.latvec.e32 << " " << GlobalC::ucell.latvec.e33 + << ucell.latvec.e31 << " " << ucell.latvec.e32 << " " << ucell.latvec.e33 << std::endl; - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < ucell.ntype; ++it) { - ofs << GlobalC::ucell.atoms[it].label << "\t"; + ofs << ucell.atoms[it].label << "\t"; } ofs << std::endl; - for (int it = 0; it < GlobalC::ucell.ntype; ++it) + for (int it = 0; it < ucell.ntype; ++it) { - ofs << GlobalC::ucell.atoms[it].na << "\t"; + ofs << ucell.atoms[it].na << "\t"; } ofs << std::endl; ofs << "Direct" << std::endl; - 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) { - ofs << GlobalC::ucell.atoms[it].taud[ia].x << " " << GlobalC::ucell.atoms[it].taud[ia].y << " " - << GlobalC::ucell.atoms[it].taud[ia].z << std::endl; + ofs << ucell.atoms[it].taud[ia].x << " " << ucell.atoms[it].taud[ia].y << " " + << ucell.atoms[it].taud[ia].z << std::endl; } } ofs << std::endl; diff --git a/source/module_io/write_wfc_r.h b/source/module_io/write_wfc_r.h index df95625bed..d2b5f76edb 100644 --- a/source/module_io/write_wfc_r.h +++ b/source/module_io/write_wfc_r.h @@ -25,7 +25,8 @@ namespace ModuleIO // write ||wfc_r|| for all k-points and all bands // Input: wfc_g[ik](ib,ig) // loop order is for(z){for(y){for(x)}} -void write_psi_r_1(const psi::Psi>& wfc_g, +void write_psi_r_1(const UnitCell& ucell, + const psi::Psi>& wfc_g, const ModulePW::PW_Basis_K* wfcpw, const std::string& folder_name, const bool& square, @@ -39,14 +40,14 @@ std::vector> cal_wfc_r(const ModulePW::PW_Basis_K* wfcpw, const int ib); // Input: chg_r[ir] -#ifdef __MPI -void write_chg_r_1(const ModulePW::PW_Basis_K* wfcpw, +void write_chg_r_1(const UnitCell& ucell, + const ModulePW::PW_Basis_K* wfcpw, const std::vector& chg_r, - const std::string& file_name, - MPI_Request& mpi_request); -#else -void write_chg_r_1(const ModulePW::PW_Basis_K* wfcpw, const std::vector& chg_r, const std::string& file_name); -#endif + const std::string& file_name + #ifdef __MPI + ,MPI_Request& mpi_request + #endif + ); } #endif