@@ -31,39 +31,38 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c
3131 ModuleBase::TITLE (" Exx_LRI" ," init" );
3232 ModuleBase::timer::tick (" Exx_LRI" , " init" );
3333
34- // if(GlobalC::exx_info.info_global.separate_loop)
35- // {
36- // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::No;
37- // Hexx_para.mixing_beta = 0;
38- // }
39- // else
40- // {
41- // if("plain"==GlobalC::CHR.mixing_mode)
42- // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Plain;
43- // else if("pulay"==GlobalC::CHR.mixing_mode)
44- // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Pulay;
45- // else
46- // throw std::invalid_argument("exx mixing error. exx_separate_loop==false, mixing_mode!=plain or pulay");
47- // Hexx_para.mixing_beta = GlobalC::CHR.mixing_beta;
48- // }
49-
50- this ->mpi_comm = mpi_comm_in;
51- this ->p_kv = &kv_in;
52- this ->orb_cutoff_ = orb.cutoffs ();
34+ // if(GlobalC::exx_info.info_global.separate_loop)
35+ // {
36+ // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::No;
37+ // Hexx_para.mixing_beta = 0;
38+ // }
39+ // else
40+ // {
41+ // if("plain"==GlobalC::CHR.mixing_mode)
42+ // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Plain;
43+ // else if("pulay"==GlobalC::CHR.mixing_mode)
44+ // Hexx_para.mixing_mode = Exx_Abfs::Parallel::Communicate::Hexx::Mixing_Mode::Pulay;
45+ // else
46+ // throw std::invalid_argument("exx mixing error. exx_separate_loop==false, mixing_mode!=plain or pulay");
47+ // Hexx_para.mixing_beta = GlobalC::CHR.mixing_beta;
48+ // }
49+
50+ this ->mpi_comm = mpi_comm_in;
51+ this ->p_kv = &kv_in;
52+ this ->orb_cutoff_ = orb.cutoffs ();
5353
5454 this ->lcaos = Exx_Abfs::Construct_Orbs::change_orbs ( orb, this ->info .kmesh_times );
5555
56- // #ifdef __MPI
57- // Exx_Abfs::Util::bcast( this->info.files_abfs, 0, this->mpi_comm );
58- // #endif
56+ // #ifdef __MPI
57+ // Exx_Abfs::Util::bcast( this->info.files_abfs, 0, this->mpi_comm );
58+ // #endif
5959
6060 const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
6161 abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom ( orb, this ->lcaos , this ->info .kmesh_times , this ->info .pca_threshold );
62- if (this ->info .files_abfs .empty ()) {
63- this ->abfs = abfs_same_atom;
64- } else {
65- this ->abfs = Exx_Abfs::IO::construct_abfs ( abfs_same_atom, orb, this ->info .files_abfs , this ->info .kmesh_times );
66- }
62+ if (this ->info .files_abfs .empty ())
63+ { this ->abfs = abfs_same_atom;}
64+ else
65+ { this ->abfs = Exx_Abfs::IO::construct_abfs ( abfs_same_atom, orb, this ->info .files_abfs , this ->info .kmesh_times ); }
6766 Exx_Abfs::Construct_Orbs::print_orbs_size (this ->abfs , GlobalV::ofs_running);
6867
6968 auto get_ccp_parameter = [this ]() -> std::map<std::string,double >
@@ -85,15 +84,14 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, c
8584 throw std::domain_error (std::string (__FILE__)+" line " +std::to_string (__LINE__)); break ;
8685 }
8786 };
88- this ->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp (this ->abfs , this ->info .ccp_type , get_ccp_parameter (), this ->info .ccp_rmesh_times );
87+ this ->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp (this ->abfs , this ->info .ccp_type , get_ccp_parameter (), this ->info .ccp_rmesh_times );
8988
9089
91- for ( size_t T=0 ; T!=this ->abfs .size (); ++T ) {
92- GlobalC::exx_info.info_ri .abfs_Lmax = std::max ( GlobalC::exx_info.info_ri .abfs_Lmax , static_cast <int >(this ->abfs [T].size ())-1 );
93- }
90+ for ( size_t T=0 ; T!=this ->abfs .size (); ++T )
91+ { GlobalC::exx_info.info_ri .abfs_Lmax = std::max ( GlobalC::exx_info.info_ri .abfs_Lmax , static_cast <int >(this ->abfs [T].size ())-1 ); }
9492
9593 this ->cv .set_orbitals (
96- orb,
94+ orb,
9795 this ->lcaos , this ->abfs , this ->abfs_ccp ,
9896 this ->info .kmesh_times , this ->info .ccp_rmesh_times );
9997
@@ -106,19 +104,17 @@ void Exx_LRI<Tdata>::cal_exx_ions(const bool write_cv)
106104 ModuleBase::TITLE (" Exx_LRI" ," cal_exx_ions" );
107105 ModuleBase::timer::tick (" Exx_LRI" , " cal_exx_ions" );
108106
109- // init_radial_table_ions( cal_atom_centres_core(atom_pairs_core_origin), atom_pairs_core_origin );
107+ // init_radial_table_ions( cal_atom_centres_core(atom_pairs_core_origin), atom_pairs_core_origin );
110108
111- // this->m_abfsabfs.init_radial_table(Rradial);
112- // this->m_abfslcaos_lcaos.init_radial_table(Rradial);
109+ // this->m_abfsabfs.init_radial_table(Rradial);
110+ // this->m_abfslcaos_lcaos.init_radial_table(Rradial);
113111
114112 std::vector<TA> atoms (GlobalC::ucell.nat );
115- for (int iat=0 ; iat<GlobalC::ucell.nat ; ++iat) {
116- atoms[iat] = iat;
117- }
113+ for (int iat=0 ; iat<GlobalC::ucell.nat ; ++iat)
114+ { atoms[iat] = iat; }
118115 std::map<TA,TatomR> atoms_pos;
119- for (int iat=0 ; iat<GlobalC::ucell.nat ; ++iat) {
120- atoms_pos[iat] = RI_Util::Vector3_to_array3 ( GlobalC::ucell.atoms [ GlobalC::ucell.iat2it [iat] ].tau [ GlobalC::ucell.iat2ia [iat] ] );
121- }
116+ for (int iat=0 ; iat<GlobalC::ucell.nat ; ++iat)
117+ { atoms_pos[iat] = RI_Util::Vector3_to_array3 ( GlobalC::ucell.atoms [ GlobalC::ucell.iat2it [iat] ].tau [ GlobalC::ucell.iat2ia [iat] ] ); }
122118 const std::array<TatomR,Ndim> latvec
123119 = {RI_Util::Vector3_to_array3 (GlobalC::ucell.a1 ),
124120 RI_Util::Vector3_to_array3 (GlobalC::ucell.a2 ),
@@ -137,7 +133,8 @@ void Exx_LRI<Tdata>::cal_exx_ions(const bool write_cv)
137133 list_As_Vs.first , list_As_Vs.second [0 ],
138134 {{" writable_Vws" ,true }});
139135 this ->cv .Vws = LRI_CV_Tools::get_CVws (Vs);
140- if (write_cv && GlobalV::MY_RANK == 0 ) { LRI_CV_Tools::write_Vs_abf (Vs, PARAM.globalv .global_out_dir + " Vs" ); }
136+ if (write_cv && GlobalV::MY_RANK == 0 )
137+ { LRI_CV_Tools::write_Vs_abf (Vs, PARAM.globalv .global_out_dir + " Vs" ); }
141138 this ->exx_lri .set_Vs (std::move (Vs), this ->info .V_threshold );
142139
143140 if (PARAM.inp .cal_force || PARAM.inp .cal_stress )
@@ -166,7 +163,8 @@ void Exx_LRI<Tdata>::cal_exx_ions(const bool write_cv)
166163 {" writable_Cws" ,true }, {" writable_dCws" ,true }, {" writable_Vws" ,false }, {" writable_dVws" ,false }});
167164 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs = std::get<0 >(Cs_dCs);
168165 this ->cv .Cws = LRI_CV_Tools::get_CVws (Cs);
169- if (write_cv && GlobalV::MY_RANK == 0 ) { LRI_CV_Tools::write_Cs_ao (Cs, PARAM.globalv .global_out_dir + " Cs" ); }
166+ if (write_cv && GlobalV::MY_RANK == 0 )
167+ { LRI_CV_Tools::write_Cs_ao (Cs, PARAM.globalv .global_out_dir + " Cs" ); }
170168 this ->exx_lri .set_Cs (std::move (Cs), this ->info .C_threshold );
171169
172170 if (PARAM.inp .cal_force || PARAM.inp .cal_stress )
@@ -185,44 +183,48 @@ void Exx_LRI<Tdata>::cal_exx_ions(const bool write_cv)
185183
186184template <typename Tdata>
187185void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
188- const Parallel_Orbitals& pv,
189- const ModuleSymmetry::Symmetry_rotation* p_symrot)
186+ const Parallel_Orbitals& pv,
187+ const ModuleSymmetry::Symmetry_rotation* p_symrot)
190188{
191189 ModuleBase::TITLE (" Exx_LRI" ," cal_exx_elec" );
192190 ModuleBase::timer::tick (" Exx_LRI" , " cal_exx_elec" );
193191
194192 const std::vector<std::tuple<std::set<TA>, std::set<TA>>> judge = RI_2D_Comm::get_2D_judge (pv);
195193
194+ if (p_symrot)
195+ { this ->exx_lri .set_symmetry (true , p_symrot->get_irreducible_sector ()); }
196+ else
197+ { this ->exx_lri .set_symmetry (false , {}); }
198+
196199 this ->Hexxs .resize (PARAM.inp .nspin );
197200 this ->Eexx = 0 ;
198- (p_symrot) ? this ->exx_lri .set_symmetry (true , p_symrot->get_irreducible_sector ()) : this ->exx_lri .set_symmetry (false , {});
199201 for (int is=0 ; is<PARAM.inp .nspin ; ++is)
200202 {
201- std::string suffix = ((PARAM.inp .cal_force || PARAM.inp .cal_stress ) ? std::to_string (is) : " " );
202-
203- this ->exx_lri .set_Ds (Ds[is], this ->info .dm_threshold , suffix);
204- this ->exx_lri .cal_Hs ({ " " ," " ,suffix });
205-
206- if (!p_symrot)
207- {
208- this ->Hexxs [is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first (
209- this ->mpi_comm , std::move (this ->exx_lri .Hs ), std::get<0 >(judge[is]), std::get<1 >(judge[is]));
210- }
211- else
212- {
213- // reduce but not repeat
214- auto Hs_a2D = this ->exx_lri .post_2D .set_tensors_map2 (this ->exx_lri .Hs );
215- // rotate locally without repeat
216- Hs_a2D = p_symrot->restore_HR (GlobalC::ucell.symm , GlobalC::ucell.atoms , GlobalC::ucell.st , ' H' , Hs_a2D);
217- // cal energy using full Hs without repeat
218- this ->exx_lri .energy = this ->exx_lri .post_2D .cal_energy (
219- this ->exx_lri .post_2D .saves [" Ds_" + suffix],
220- this ->exx_lri .post_2D .set_tensors_map2 (Hs_a2D));
221- // get repeated full Hs for abacus
222- this ->Hexxs [is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first (
223- this ->mpi_comm , std::move (Hs_a2D), std::get<0 >(judge[is]), std::get<1 >(judge[is]));
224- }
225- this ->Eexx += std::real (this ->exx_lri .energy );
203+ const std::string suffix = ((PARAM.inp .cal_force || PARAM.inp .cal_stress ) ? std::to_string (is) : " " );
204+
205+ this ->exx_lri .set_Ds (Ds[is], this ->info .dm_threshold , suffix);
206+ this ->exx_lri .cal_Hs ({ " " ," " ,suffix });
207+
208+ if (!p_symrot)
209+ {
210+ this ->Hexxs [is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first (
211+ this ->mpi_comm , std::move (this ->exx_lri .Hs ), std::get<0 >(judge[is]), std::get<1 >(judge[is]));
212+ }
213+ else
214+ {
215+ // reduce but not repeat
216+ auto Hs_a2D = this ->exx_lri .post_2D .set_tensors_map2 (this ->exx_lri .Hs );
217+ // rotate locally without repeat
218+ Hs_a2D = p_symrot->restore_HR (GlobalC::ucell.symm , GlobalC::ucell.atoms , GlobalC::ucell.st , ' H' , Hs_a2D);
219+ // cal energy using full Hs without repeat
220+ this ->exx_lri .energy = this ->exx_lri .post_2D .cal_energy (
221+ this ->exx_lri .post_2D .saves [" Ds_" + suffix],
222+ this ->exx_lri .post_2D .set_tensors_map2 (Hs_a2D));
223+ // get repeated full Hs for abacus
224+ this ->Hexxs [is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first (
225+ this ->mpi_comm , std::move (Hs_a2D), std::get<0 >(judge[is]), std::get<1 >(judge[is]));
226+ }
227+ this ->Eexx += std::real (this ->exx_lri .energy );
226228 post_process_Hexx (this ->Hexxs [is]);
227229 }
228230 this ->Eexx = post_process_Eexx (this ->Eexx );
@@ -245,8 +247,8 @@ template<typename Tdata>
245247double Exx_LRI<Tdata>::post_process_Eexx(const double & Eexx_in) const
246248{
247249 ModuleBase::TITLE (" Exx_LRI" ," post_process_Eexx" );
248- const double SPIN_multiple = std::map<int , double >{ {1 ,2 }, {2 ,1 }, {4 ,1 } }.at (PARAM.inp .nspin ); // why?
249- const double frac = -SPIN_multiple;
250+ const double SPIN_multiple = std::map<int , double >{ {1 ,2 }, {2 ,1 }, {4 ,1 } }.at (PARAM.inp .nspin ); // why?
251+ const double frac = -SPIN_multiple;
250252 return frac * Eexx_in;
251253}
252254
@@ -280,8 +282,7 @@ void Exx_LRI<Tdata>::cal_exx_force()
280282 for (std::size_t idim=0 ; idim<Ndim; ++idim) {
281283 for (const auto &force_item : this ->exx_lri .force [idim]) {
282284 this ->force_exx (force_item.first , idim) += std::real (force_item.second );
283- }
284- }
285+ } }
285286 }
286287
287288 const double SPIN_multiple = std::map<int ,double >{{1 ,2 }, {2 ,1 }, {4 ,1 }}.at (PARAM.inp .nspin ); // why?
@@ -304,8 +305,7 @@ void Exx_LRI<Tdata>::cal_exx_stress()
304305 for (std::size_t idim0=0 ; idim0<Ndim; ++idim0) {
305306 for (std::size_t idim1=0 ; idim1<Ndim; ++idim1) {
306307 this ->stress_exx (idim0,idim1) += std::real (this ->exx_lri .stress (idim0,idim1));
307- }
308- }
308+ } }
309309 }
310310
311311 const double SPIN_multiple = std::map<int ,double >{{1 ,2 }, {2 ,1 }, {4 ,1 }}.at (PARAM.inp .nspin ); // why?
@@ -318,16 +318,15 @@ void Exx_LRI<Tdata>::cal_exx_stress()
318318template <typename Tdata>
319319std::vector<std::vector<int >> Exx_LRI<Tdata>::get_abfs_nchis() const
320320{
321- std::vector<std::vector<int >> abfs_nchis;
322- for (const auto & abfs_T : this ->abfs )
323- {
324- std::vector<int > abfs_nchi_T;
325- for (const auto & abfs_L : abfs_T) {
326- abfs_nchi_T.push_back (abfs_L.size ());
327- }
328- abfs_nchis.push_back (abfs_nchi_T);
329- }
330- return abfs_nchis;
321+ std::vector<std::vector<int >> abfs_nchis;
322+ for (const auto & abfs_T : this ->abfs )
323+ {
324+ std::vector<int > abfs_nchi_T;
325+ for (const auto & abfs_L : abfs_T)
326+ { abfs_nchi_T.push_back (abfs_L.size ()); }
327+ abfs_nchis.push_back (abfs_nchi_T);
328+ }
329+ return abfs_nchis;
331330}
332331
333332#endif
0 commit comments