14 using TAC = std::pair<int, std::array<int, 3>>;
17 template <
typename Tdata,
typename TR>
18 void reallocate_hcontainer(
const std::vector<std::map<
int, std::map<
TAC, RI::Tensor<Tdata>>>>& Hexxs,
20 const RI::Cell_Nearest<int, int, 3, double, 3>*
const cell_nearest)
22 auto* pv = hR->get_paraV();
23 bool need_allocate =
false;
24 for (
auto& Htmp1 : Hexxs[0])
26 const int& iat0 = Htmp1.first;
27 for (
auto& Htmp2 : Htmp1.second)
29 const int& iat1 = Htmp2.first.first;
30 if (pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0)
34 cell_nearest->get_cell_nearest_discrete(iat0, iat1, Htmp2.first.second)
35 : Htmp2.first.second));
36 BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.
x, R.
y, R.
z);
40 AtomPair<TR> tmp(iat0, iat1, R.
x, R.
y, R.
z, pv);
46 if (need_allocate) { hR->allocate(
nullptr,
true); }
50 template <
typename TR>
51 void reallocate_hcontainer(
const int nat, HContainer<TR>* hR,
52 const std::array<int, 3>& Rs_period,
53 const RI::Cell_Nearest<int, int, 3, double, 3>*
const cell_nearest)
55 auto* pv = hR->get_paraV();
57 bool need_allocate =
false;
58 for (
int iat0 = 0;iat0 < nat;++iat0)
60 for (
int iat1 = 0;iat1 < nat;++iat1)
64 if(pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0)
70 cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell)
72 BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.
x, R.
y, R.
z);
77 AtomPair<TR> tmp(iat0, iat1, R.
x, R.
y, R.
z, pv);
84 if (need_allocate) { hR->allocate(
nullptr,
true);}
87template <
typename TK,
typename TR>
88OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
92 std::vector<std::map<
int, std::map<
TAC, RI::Tensor<double>>>>* Hexxd_in,
93 std::vector<std::map<
int, std::map<
TAC, RI::Tensor<std::complex<double>>>>>* Hexxc_in,
94 Add_Hexx_Type add_hexx_type_in,
96 int* two_level_step_in,
97 const bool restart_in)
98 : OperatorLCAO<TK, TR>(hsk_in, kv_in.kvec_d, hR_in),
103 add_hexx_type(add_hexx_type_in),
105 two_level_step(two_level_step_in),
109 this->cal_type = calculation_type::lcao_exx;
114 auto file_name_list_csr = []() -> std::vector<std::string>
116 std::vector<std::string> file_name_list;
121 return file_name_list;
123 auto file_name_list_cereal = []() -> std::vector<std::string>
125 std::vector<std::string> file_name_list;
127 { file_name_list.push_back(
"HexxR_" + std::to_string(irank) ); }
128 return file_name_list;
130 auto check_exist = [](
const std::vector<std::string> &file_name_list) ->
bool
132 for (
const std::string &file_name : file_name_list)
134 std::ifstream ifs(file_name);
141 std::cout<<
" Attention: The number of MPI processes must be strictly identical between SCF and NSCF when computing exact-exchange."<<std::endl;
142 if (check_exist(file_name_list_csr()))
149 if (this->add_hexx_type == Add_Hexx_Type::R)
150 { reallocate_hcontainer(*Hexxd, this->hR); }
155 if (this->add_hexx_type == Add_Hexx_Type::R)
156 { reallocate_hcontainer(*Hexxc, this->hR); }
159 else if (check_exist(file_name_list_cereal()))
163 std::ifstream ifs(file_name_exx_cereal, std::ios::binary);
169 if (this->add_hexx_type == Add_Hexx_Type::R)
170 { reallocate_hcontainer(*Hexxd, this->hR); }
175 if (this->add_hexx_type == Add_Hexx_Type::R)
176 { reallocate_hcontainer(*Hexxc, this->hR); }
183 this->use_cell_nearest =
false;
187 if (this->add_hexx_type == Add_Hexx_Type::R)
191 std::fmod(this->kv.get_koffset(1), 1.0), std::fmod(this->kv.get_koffset(2), 1.0)).
norm() < 1e-10);
193 const std::array<int, 3> Rs_period = { this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2] };
194 if (this->use_cell_nearest)
196 this->cell_nearest = init_cell_nearest(ucell, Rs_period);
197 reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest);
199 else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); }
205 assert(this->two_level_step !=
nullptr);
207 if (this->add_hexx_type == Add_Hexx_Type::k)
210 if (std::is_same<TK, double>::value)
212 this->Hexxd_k_load.resize(this->kv.get_nks());
213 for (
int ik = 0; ik < this->kv.get_nks(); ik++)
219 if (!this->restart) {
break; }
224 this->Hexxc_k_load.resize(this->kv.get_nks());
225 for (
int ik = 0; ik < this->kv.get_nks(); ik++)
231 if (!this->restart) {
break; }
235 else if (this->add_hexx_type == Add_Hexx_Type::R)
242 std::ifstream ifs(restart_HR_path +
"_" + std::to_string(is) +
".csr");
243 if (!ifs) { all_exist = 0;
break; }
248 MPI_Allreduce(MPI_IN_PLACE, &all_exist, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
264 std::ifstream ifs(restart_HR_path_cereal, std::ios::binary);
265 int all_exist_cereal = ifs ? 1 : 0;
267 MPI_Allreduce(MPI_IN_PLACE, &all_exist_cereal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
269 if (!all_exist_cereal)
286 if (!this->restart) {
287 std::cout <<
"WARNING: Hexx not found, restart from the non-exx loop." << std::endl
288 <<
"If the loaded charge density is EXX-solved, this may lead to poor convergence." << std::endl;
295template<
typename TK,
typename TR>
296void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
302 && this->two_level_step !=
nullptr && *this->two_level_step == 0
308 if (this->add_hexx_type == Add_Hexx_Type::k) {
return; }
319 *this->hR->get_paraV(),
322 this->use_cell_nearest ? &this->cell_nearest : nullptr);
330 *this->hR->get_paraV(),
333 this->use_cell_nearest ? &this->cell_nearest : nullptr);
336 if (
PARAM.
inp.
nspin == 2) { this->current_spin = 1 - this->current_spin; }
339template<
typename TK,
typename TR>
340void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(
int ik)
344 if (
PARAM.
inp.
calculation !=
"nscf" && this->two_level_step !=
nullptr && *this->two_level_step == 0 && !this->restart) {
return; }
346 if (this->add_hexx_type == Add_Hexx_Type::R) {
throw std::invalid_argument(
"Set Add_Hexx_Type::k sto call OperatorEXX::contributeHk()."); }
350 if (this->restart && this->two_level_step !=
nullptr)
352 if (*this->two_level_step == 0)
354 this->add_loaded_Hexx(ik);
359 if (this->Hexxd_k_load.size() > 0)
361 this->Hexxd_k_load.clear();
362 this->Hexxd_k_load.shrink_to_fit();
364 else if (this->Hexxc_k_load.size() > 0)
366 this->Hexxc_k_load.clear();
367 this->Hexxc_k_load.shrink_to_fit();
380 *this->hR->get_paraV(),
381 this->hsk->get_hk());
389 *this->hR->get_paraV(),
390 this->hsk->get_hk());
Definition abfs-vector3_order.h:16
3 elements vector
Definition vector3.h:22
T x
Definition vector3.h:24
T norm(void) const
Get the norm of a Vector3.
Definition vector3.h:187
T y
Definition vector3.h:25
T z
Definition vector3.h:26
int64_t get_local_size() const
number of local matrix elements
Definition parallel_2d.h:39
Definition parallel_orbitals.h:9
const Input_para & inp
Definition parameter.h:26
const System_para & globalv
Definition parameter.h:30
std::string folder
Definition restart.h:31
Info_Load info_load
Definition restart.h:29
bool load_disk(const std::string label, const int index, const int size, T *data, const bool error_quit=true) const
Definition restart.h:43
static int get_func_type()
Definition xc_functional.h:67
Exx_Info exx_info
Definition test_xc.cpp:29
Restart restart
Definition for_testing_input_conv.h:255
int istep
Definition ions_move_basic.cpp:12
void WARNING_QUIT(const std::string &, const std::string &)
Combine the functions of WARNING and QUIT.
Definition test_delley.cpp:14
void TITLE(const std::string &class_name, const std::string &function_name, const bool disable)
Definition tool_title.cpp:18
void read_Hexxs_cereal(const std::string &file_name, std::vector< std::map< int, std::map< TAC, RI::Tensor< Tdata > > > > &Hexxs)
read Hexxs in cereal format
Definition restart_exx_csr.hpp:64
void read_Hexxs_csr(const std::string &file_name, const UnitCell &ucell, const int nspin, const int nbasis, std::vector< std::map< int, std::map< TAC, RI::Tensor< Tdata > > > > &Hexxs)
read Hexxs in CSR format
Definition restart_exx_csr.hpp:13
void add_HexxR(const int current_spin, const double alpha, const std::vector< std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > > &Hs, const Parallel_Orbitals &pv, const int npol, hamilt::HContainer< TR > &HlocR, const RI::Cell_Nearest< int, int, 3, double, 3 > *const cell_nearest=nullptr)
Definition RI_2D_Comm.hpp:238
void add_Hexx(const UnitCell &ucell, const K_Vectors &kv, const int ik, const double alpha, const std::vector< std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > > &Hs, const Parallel_Orbitals &pv, TK *hk)
Definition RI_2D_Comm.hpp:125
std::vector< std::array< Tcell, Ndim > > get_Born_von_Karmen_cells(const std::array< Tcell, Ndim > &Born_von_Karman_period)
Definition RI_Util.hpp:34
ModuleBase::Vector3< Tcell > array3_to_Vector3(const std::array< Tcell, 3 > &v)
Definition RI_Util.h:38
Parameter PARAM
Definition parameter.cpp:3
std::pair< int, TC > TAC
Definition ri_cv_io_test.cpp:10
bool cal_exx
Definition exx_info.h:15
Exx_Info_Global info_global
Definition exx_info.h:38
bool load_H_finish
Definition restart.h:26
std::string global_readin_dir
global readin directory
Definition system_parameter.h:43
int npol
number of polarization
Definition system_parameter.h:52
int nlocal
total number of local basis.
Definition system_parameter.h:23
int myrank
Definition system_parameter.h:11
int nproc
Definition system_parameter.h:12