ABACUS develop
Atomic-orbital Based Ab-initio Computation at UStc
Loading...
Searching...
No Matches
op_exx_lcao.hpp
Go to the documentation of this file.
1#ifndef OPEXXLCAO_HPP
2#define OPEXXLCAO_HPP
3#ifdef __EXX
4
5#include "op_exx_lcao.h"
11
12namespace hamilt
13{
14 using TAC = std::pair<int, std::array<int, 3>>;
15
16 // allocate according to the read-in HexxR, used in nscf
17 template <typename Tdata, typename TR>
18 void reallocate_hcontainer(const std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& Hexxs,
19 HContainer<TR>* hR,
20 const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest)
21 {
22 auto* pv = hR->get_paraV();
23 bool need_allocate = false;
24 for (auto& Htmp1 : Hexxs[0])
25 {
26 const int& iat0 = Htmp1.first;
27 for (auto& Htmp2 : Htmp1.second)
28 {
29 const int& iat1 = Htmp2.first.first;
30 if (pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0)
31 {
33 (cell_nearest ?
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);
37 if (HlocR == nullptr)
38 { // add R to HContainer
39 need_allocate = true;
40 AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, pv);
41 hR->insert_pair(tmp);
42 }
43 }
44 }
45 }
46 if (need_allocate) { hR->allocate(nullptr, true); }
47 }
48
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)
54 {
55 auto* pv = hR->get_paraV();
56 auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period);
57 bool need_allocate = false;
58 for (int iat0 = 0;iat0 < nat;++iat0)
59 {
60 for (int iat1 = 0;iat1 < nat;++iat1)
61 {
62 // complete the atom pairs that has orbitals in this processor but not in hR due to the adj_list
63 // but adj_list is not enought for EXX, which is more nonlocal than Nonlocal
64 if(pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0)
65 {
66 for (auto& cell : Rs)
67 {
69 (cell_nearest ?
70 cell_nearest->get_cell_nearest_discrete(iat0, iat1, cell)
71 : cell));
72 BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
73
74 if (HlocR == nullptr)
75 { // add R to HContainer
76 need_allocate = true;
77 AtomPair<TR> tmp(iat0, iat1, R.x, R.y, R.z, pv);
78 hR->insert_pair(tmp);
79 }
80 }
81 }
82 }
83 }
84 if (need_allocate) { hR->allocate(nullptr, true);}
85 }
86
87template <typename TK, typename TR>
88OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
89 HContainer<TR>*hR_in,
90 const UnitCell& ucell_in,
91 const K_Vectors& kv_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,
95 const int istep,
96 int* two_level_step_in,
97 const bool restart_in)
98 : OperatorLCAO<TK, TR>(hsk_in, kv_in.kvec_d, hR_in),
99 ucell(ucell_in),
100 kv(kv_in),
101 Hexxd(Hexxd_in),
102 Hexxc(Hexxc_in),
103 add_hexx_type(add_hexx_type_in),
104 istep(istep),
105 two_level_step(two_level_step_in),
106 restart(restart_in)
107{
108 ModuleBase::TITLE("OperatorEXX", "OperatorEXX");
109 this->cal_type = calculation_type::lcao_exx;
110 const Parallel_Orbitals* const pv = hR_in->get_paraV();
111
113 { // if nscf, read HexxR first and reallocate hR according to the read-in HexxR
114 auto file_name_list_csr = []() -> std::vector<std::string>
115 {
116 std::vector<std::string> file_name_list;
117 for (int irank=0; irank<PARAM.globalv.nproc; ++irank) {
118 for (int is=0;is<PARAM.inp.nspin;++is) {
119 file_name_list.push_back( PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(irank) + "_" + std::to_string(is) + ".csr" );
120 } }
121 return file_name_list;
122 };
123 auto file_name_list_cereal = []() -> std::vector<std::string>
124 {
125 std::vector<std::string> file_name_list;
126 for (int irank=0; irank<PARAM.globalv.nproc; ++irank)
127 { file_name_list.push_back( "HexxR_" + std::to_string(irank) ); }
128 return file_name_list;
129 };
130 auto check_exist = [](const std::vector<std::string> &file_name_list) -> bool
131 {
132 for (const std::string &file_name : file_name_list)
133 {
134 std::ifstream ifs(file_name);
135 if (!ifs.is_open())
136 { return false; }
137 }
138 return true;
139 };
140
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()))
143 {
144 const std::string file_name_exx_csr = PARAM.globalv.global_readin_dir + "HexxR" + std::to_string(PARAM.globalv.myrank);
145 // Read HexxR in CSR format
146 if (GlobalC::exx_info.info_ri.real_number)
147 {
148 ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
149 if (this->add_hexx_type == Add_Hexx_Type::R)
150 { reallocate_hcontainer(*Hexxd, this->hR); }
151 }
152 else
153 {
154 ModuleIO::read_Hexxs_csr(file_name_exx_csr, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
155 if (this->add_hexx_type == Add_Hexx_Type::R)
156 { reallocate_hcontainer(*Hexxc, this->hR); }
157 }
158 }
159 else if (check_exist(file_name_list_cereal()))
160 {
161 // Read HexxR in binary format (old version)
162 const std::string file_name_exx_cereal = PARAM.globalv.global_readin_dir + "HexxR_" + std::to_string(PARAM.globalv.myrank);
163 std::ifstream ifs(file_name_exx_cereal, std::ios::binary);
164 if (!ifs)
165 { ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file < " + file_name_exx_cereal + " >."); }
166 if (GlobalC::exx_info.info_ri.real_number)
167 {
168 ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxd);
169 if (this->add_hexx_type == Add_Hexx_Type::R)
170 { reallocate_hcontainer(*Hexxd, this->hR); }
171 }
172 else
173 {
174 ModuleIO::read_Hexxs_cereal(file_name_exx_cereal, *Hexxc);
175 if (this->add_hexx_type == Add_Hexx_Type::R)
176 { reallocate_hcontainer(*Hexxc, this->hR); }
177 }
178 }
179 else
180 {
181 ModuleBase::WARNING_QUIT("OperatorEXX", "Can't open EXX file in " + PARAM.globalv.global_readin_dir);
182 }
183 this->use_cell_nearest = false;
184 }
185 else
186 { // if scf and Add_Hexx_Type::R, init cell_nearest and reallocate hR according to BvK cells
187 if (this->add_hexx_type == Add_Hexx_Type::R)
188 {
189 // if k points has no shift, use cell_nearest to reduce the memory cost
190 this->use_cell_nearest = (ModuleBase::Vector3<double>(std::fmod(this->kv.get_koffset(0), 1.0),
191 std::fmod(this->kv.get_koffset(1), 1.0), std::fmod(this->kv.get_koffset(2), 1.0)).norm() < 1e-10);
192
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)
195 {
196 this->cell_nearest = init_cell_nearest(ucell, Rs_period);
197 reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest);
198 }
199 else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); }
200 }
201
202 if (this->restart)
203 {
205 assert(this->two_level_step != nullptr);
206
207 if (this->add_hexx_type == Add_Hexx_Type::k)
208 {
210 if (std::is_same<TK, double>::value)
211 {
212 this->Hexxd_k_load.resize(this->kv.get_nks());
213 for (int ik = 0; ik < this->kv.get_nks(); ik++)
214 {
215 this->Hexxd_k_load[ik].resize(pv->get_local_size(), 0.0);
217 "Hexx", ik,
218 pv->get_local_size(), this->Hexxd_k_load[ik].data(), false);
219 if (!this->restart) { break; }
220 }
221 }
222 else
223 {
224 this->Hexxc_k_load.resize(this->kv.get_nks());
225 for (int ik = 0; ik < this->kv.get_nks(); ik++)
226 {
227 this->Hexxc_k_load[ik].resize(pv->get_local_size(), 0.0);
229 "Hexx", ik,
230 pv->get_local_size(), this->Hexxc_k_load[ik].data(), false);
231 if (!this->restart) { break; }
232 }
233 }
234 }
235 else if (this->add_hexx_type == Add_Hexx_Type::R)
236 {
237 // read in Hexx(R)
238 const std::string restart_HR_path = GlobalC::restart.folder + "HexxR" + std::to_string(PARAM.globalv.myrank);
239 int all_exist = 1;
240 for (int is = 0; is < PARAM.inp.nspin; ++is)
241 {
242 std::ifstream ifs(restart_HR_path + "_" + std::to_string(is) + ".csr");
243 if (!ifs) { all_exist = 0; break; }
244 }
245// Add MPI communication to synchronize all_exist across processes
246#ifdef __MPI
247 // don't read in any files if one of the processes doesn't have it
248 MPI_Allreduce(MPI_IN_PLACE, &all_exist, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
249#endif
250 if (all_exist)
251 {
252 // Read HexxR in CSR format
253 if (GlobalC::exx_info.info_ri.real_number) {
254 ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxd);
255 }
256 else {
257 ModuleIO::read_Hexxs_csr(restart_HR_path, ucell, PARAM.inp.nspin, PARAM.globalv.nlocal, *Hexxc);
258 }
259 }
260 else
261 {
262 // Read HexxR in binary format (old version)
263 const std::string restart_HR_path_cereal = GlobalC::restart.folder + "HexxR_" + std::to_string(PARAM.globalv.myrank);
264 std::ifstream ifs(restart_HR_path_cereal, std::ios::binary);
265 int all_exist_cereal = ifs ? 1 : 0;
266#ifdef __MPI
267 MPI_Allreduce(MPI_IN_PLACE, &all_exist_cereal, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
268#endif
269 if (!all_exist_cereal)
270 {
271 //no HexxR file in CSR or binary format
272 this->restart = false;
273 }
274 else
275 {
276 if (GlobalC::exx_info.info_ri.real_number) {
277 ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxd);
278 }
279 else {
280 ModuleIO::read_Hexxs_cereal(restart_HR_path_cereal, *Hexxc);
281 }
282 }
283 }
284 }
285
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;
289 }
291 }
292 }
293}
294
295template<typename TK, typename TR>
296void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHR()
297{
298 ModuleBase::TITLE("OperatorEXX", "contributeHR");
299 // Peize Lin add 2016-12-03
300 if (this->istep == 0
301 && PARAM.inp.calculation != "nscf"
302 && this->two_level_step != nullptr && *this->two_level_step == 0
303 && PARAM.inp.init_wfc != "file"
304 && !this->restart)
305 {
306 return;
307 } //in the non-exx loop, do nothing
308 if (this->add_hexx_type == Add_Hexx_Type::k) { return; }
309
311 {
312 // add H(R) normally
313 if (GlobalC::exx_info.info_ri.real_number)
314 {
316 this->current_spin,
317 GlobalC::exx_info.info_global.hybrid_alpha,
318 *this->Hexxd,
319 *this->hR->get_paraV(),
321 *this->hR,
322 this->use_cell_nearest ? &this->cell_nearest : nullptr);
323 }
324 else
325 {
327 this->current_spin,
328 GlobalC::exx_info.info_global.hybrid_alpha,
329 *this->Hexxc,
330 *this->hR->get_paraV(),
332 *this->hR,
333 this->use_cell_nearest ? &this->cell_nearest : nullptr);
334 }
335 }
336 if (PARAM.inp.nspin == 2) { this->current_spin = 1 - this->current_spin; }
337}
338
339template<typename TK, typename TR>
340void OperatorEXX<OperatorLCAO<TK, TR>>::contributeHk(int ik)
341{
342 ModuleBase::TITLE("OperatorEXX", "constributeHR");
343 // Peize Lin add 2016-12-03
344 if (PARAM.inp.calculation != "nscf" && this->two_level_step != nullptr && *this->two_level_step == 0 && !this->restart) { return; } //in the non-exx loop, do nothing
345
346 if (this->add_hexx_type == Add_Hexx_Type::R) { throw std::invalid_argument("Set Add_Hexx_Type::k sto call OperatorEXX::contributeHk()."); }
347
349 {
350 if (this->restart && this->two_level_step != nullptr)
351 {
352 if (*this->two_level_step == 0)
353 {
354 this->add_loaded_Hexx(ik);
355 return;
356 }
357 else // clear loaded Hexx and release memory
358 {
359 if (this->Hexxd_k_load.size() > 0)
360 {
361 this->Hexxd_k_load.clear();
362 this->Hexxd_k_load.shrink_to_fit();
363 }
364 else if (this->Hexxc_k_load.size() > 0)
365 {
366 this->Hexxc_k_load.clear();
367 this->Hexxc_k_load.shrink_to_fit();
368 }
369 }
370 }
371 // cal H(k) from H(R) normally
372
373 if (GlobalC::exx_info.info_ri.real_number) {
375 ucell,
376 this->kv,
377 ik,
378 GlobalC::exx_info.info_global.hybrid_alpha,
379 *this->Hexxd,
380 *this->hR->get_paraV(),
381 this->hsk->get_hk());
382 } else {
384 ucell,
385 this->kv,
386 ik,
387 GlobalC::exx_info.info_global.hybrid_alpha,
388 *this->Hexxc,
389 *this->hR->get_paraV(),
390 this->hsk->get_hk());
391}
392 }
393}
394
395} // namespace hamilt
396#endif // __EXX
397#endif // OPEXXLCAO_HPP
Definition abfs-vector3_order.h:16
Definition klist.h:13
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
Definition unitcell.h:16
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
Definition hamilt.h:12
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
std::string calculation
Definition input_parameter.h:19
std::string init_wfc
"file","atomic","random"
Definition input_parameter.h:47
int nspin
LDA ; LSDA ; non-linear spin.
Definition input_parameter.h:84
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