ABACUS develop
Atomic-orbital Based Ab-initio Computation at UStc
Loading...
Searching...
No Matches
Exx_LRI.hpp
Go to the documentation of this file.
1//=======================
2// AUTHOR : Peize Lin
3#include "source_io/module_parameter/parameter.h"
4// DATE : 2022-08-17
5//=======================
6
7#ifndef EXX_LRI_HPP
8#define EXX_LRI_HPP
9
10#include "Exx_LRI.h"
11#include "RI_2D_Comm.h"
12#include "RI_Util.h"
13#include "source_lcao/module_ri/exx_abfs-construct_orbs.h"
14#include "source_lcao/module_ri/exx_abfs-io.h"
15#include "source_lcao/module_ri/conv_coulomb_pot_k.h"
16#include "source_base/tool_title.h"
17#include "source_base/timer.h"
18#include "source_lcao/module_ri/serialization_cereal.h"
19#include "source_lcao/module_ri/Mix_DMk_2D.h"
20#include "source_basis/module_ao/parallel_orbitals.h"
21
22#include <RI/distribute/Distribute_Equally.h>
23#include <RI/global/Map_Operator-3.h>
24
25#include <fstream>
26#include <stdexcept>
27#include <string>
28
29#if defined(__GLIBC__)
30#include <malloc.h>
31#endif
32
33namespace ExxLriDetail
34{
35using CoulombParam
36 = std::map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string, std::string>>>;
37
38inline void trim_malloc_cache()
39{
40#if defined(__GLIBC__)
41 malloc_trim(0);
42#endif
43}
44
45inline double default_spencer_rcut(const UnitCell& ucell, const K_Vectors& kv)
46{
47 return std::pow(0.75 * kv.get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0);
48}
49
50inline CoulombParam build_center2_cut_coulomb_param(const CoulombParam& coulomb_param,
51 const UnitCell& ucell,
52 const K_Vectors& kv,
53 bool* synthesized_rcut = nullptr)
54{
55 CoulombParam center2_param = RI_Util::update_coulomb_param(coulomb_param, ucell, &kv);
56 const double fallback_rcut = default_spencer_rcut(ucell, kv);
57 bool used_fallback_rcut = false;
58
59 for (auto& param_list: center2_param)
60 {
61 if (param_list.first != Conv_Coulomb_Pot_K::Coulomb_Type::Fock)
62 {
63 continue;
64 }
65 for (auto& param: param_list.second)
66 {
67 auto rcut_it = param.find("Rcut");
68 if (rcut_it == param.end() || rcut_it->second.empty())
69 {
70 param["Rcut"] = ModuleBase::GlobalFunc::TO_STRING(fallback_rcut);
71 used_fallback_rcut = true;
72 }
73 }
74 }
75
76 if (synthesized_rcut != nullptr)
77 {
78 *synthesized_rcut = used_fallback_rcut;
79 }
80 return center2_param;
81}
82}
83
84template<typename Tdata>
85void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
86 const UnitCell &ucell,
87 const K_Vectors &kv_in,
88 const LCAO_Orbitals& orb)
89{
90 ModuleBase::TITLE("Exx_LRI","init");
91 ModuleBase::timer::start("Exx_LRI", "init");
92
93 this->mpi_comm = mpi_comm_in;
94 this->p_kv = &kv_in;
95 this->orb_cutoff_ = orb.cutoffs();
96
97 this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times );
98 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos);
99
100 const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
101 abfs_same_atom = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell, orb, this->lcaos, this->info.kmesh_times, this->info.pca_threshold );
102 if(this->info.files_abfs.empty())
103 { this->abfs = abfs_same_atom;}
104 else
105 { this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); }
106 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs);
107 Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
108
109 for( size_t T=0; T!=this->abfs.size(); ++T )
110 { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size())-1 ); }
111
112 this->exx_objs.clear();
113 this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv);
114
115 this->MGT = std::make_shared<ORB_gaunt_table>();
116 for(const auto &settings_list : this->coulomb_settings)
117 {
118 this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, settings_list.second.second, this->info.ccp_rmesh_times);
119 this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb,
120 this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp,
121 this->info.kmesh_times, this->MGT, settings_list.second.first );
122 if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald)
123 {
124 this->exx_objs[settings_list.first].evq.init(ucell, orb,
125 this->mpi_comm, this->p_kv, this->lcaos, this->abfs,
126 settings_list.second.second, this->MGT, this->info.ccp_rmesh_times, this->info.kmesh_times);
127 }
128 }
129
130 ModuleBase::timer::end("Exx_LRI", "init");
131}
132
133template<typename Tdata>
134void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
135 const UnitCell &ucell,
136 const K_Vectors &kv_in,
137 const LCAO_Orbitals& orb,
138 const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& abfs_in)
139{
140 ModuleBase::TITLE("Exx_LRI","init");
141 ModuleBase::timer::start("Exx_LRI", "init");
142
143 this->mpi_comm = mpi_comm_in;
144 this->p_kv = &kv_in;
145 this->orb_cutoff_ = orb.cutoffs();
146
147 this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs( orb, this->info.kmesh_times );
148 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos);
149
150 this->abfs = abfs_in;
151 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs);
152 Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
153
154 for( size_t T=0; T!=this->abfs.size(); ++T )
155 { GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size())-1 ); }
156
157 this->exx_objs.clear();
158 this->coulomb_settings = RI_Util::update_coulomb_settings(this->info.coulomb_param, ucell, this->p_kv);
159
160 this->MGT = std::make_shared<ORB_gaunt_table>();
161 for(const auto &settings_list : this->coulomb_settings)
162 {
163 this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, settings_list.second.second, this->info.ccp_rmesh_times);
164 this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb,
165 this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp,
166 this->info.kmesh_times, this->MGT, settings_list.second.first );
167 if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald)
168 {
169 this->exx_objs[settings_list.first].evq.init(ucell, orb,
170 this->mpi_comm, this->p_kv, this->lcaos, this->abfs,
171 settings_list.second.second, this->MGT, this->info.ccp_rmesh_times, this->info.kmesh_times);
172 }
173 }
174
175 ModuleBase::timer::end("Exx_LRI", "init");
176}
177
178template <typename Tdata>
179void Exx_LRI<Tdata>::init_spencer(const MPI_Comm& mpi_comm_in,
180 const UnitCell& ucell,
181 const K_Vectors& kv_in,
182 const LCAO_Orbitals& orb)
183{
184 ModuleBase::TITLE("Exx_LRI", "init_spencer");
185 ModuleBase::timer::start("Exx_LRI", "init_spencer");
186
187 this->mpi_comm = mpi_comm_in;
188 this->p_kv = &kv_in;
189 this->orb_cutoff_ = orb.cutoffs();
190
191 this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs(orb, this->info.kmesh_times);
192 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos);
193
194 const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_same_atom
195 = Exx_Abfs::Construct_Orbs::abfs_same_atom(ucell,
196 orb,
197 this->lcaos,
198 this->info.kmesh_times,
199 this->info.pca_threshold);
200 if (this->info.files_abfs.empty())
201 {
202 this->abfs = abfs_same_atom;
203 }
204 else
205 {
206 this->abfs = Exx_Abfs::IO::construct_abfs(abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times);
207 }
208 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs);
209 Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
210
211 for (size_t T = 0; T != this->abfs.size(); ++T)
212 {
213 GlobalC::exx_info.info_ri.abfs_Lmax
214 = std::max(GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size()) - 1);
215 }
216
217 this->exx_objs.clear();
218 this->coulomb_settings.clear();
219 this->coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2]
220 = std::make_pair(true,
221 ExxLriDetail::build_center2_cut_coulomb_param(
222 this->info.coulomb_param, ucell, kv_in));
223
224 this->MGT = std::make_shared<ORB_gaunt_table>();
225 const auto center2_settings = this->coulomb_settings.find(Conv_Coulomb_Pot_K::Coulomb_Method::Center2);
226 if (center2_settings == this->coulomb_settings.end())
227 {
228 throw std::invalid_argument("Exx_LRI::init_spencer failed to prepare Center2 settings.");
229 }
230
231 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp_spencer(
232 this->abfs,
233 center2_settings->second.second,
234 this->info.ccp_rmesh_times);
235 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.set_orbitals(
236 ucell,
237 orb,
238 this->lcaos,
239 this->abfs,
240 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].abfs_ccp,
241 this->info.kmesh_times,
242 this->MGT,
243 center2_settings->second.first);
244
245 ModuleBase::timer::end("Exx_LRI", "init_spencer");
246}
247
248template <typename Tdata>
249void Exx_LRI<Tdata>::init_spencer(const MPI_Comm& mpi_comm_in,
250 const UnitCell& ucell,
251 const K_Vectors& kv_in,
252 const LCAO_Orbitals& orb,
253 const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& abfs_in)
254{
255 ModuleBase::TITLE("Exx_LRI", "init_spencer");
256 ModuleBase::timer::start("Exx_LRI", "init_spencer");
257
258 this->mpi_comm = mpi_comm_in;
259 this->p_kv = &kv_in;
260 this->orb_cutoff_ = orb.cutoffs();
261
262 this->lcaos = Exx_Abfs::Construct_Orbs::change_orbs(orb, this->info.kmesh_times);
263 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->lcaos);
264
265 this->abfs = abfs_in;
266 Exx_Abfs::Construct_Orbs::filter_empty_orbs(this->abfs);
267 Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
268
269 for (size_t T = 0; T != this->abfs.size(); ++T)
270 {
271 GlobalC::exx_info.info_ri.abfs_Lmax
272 = std::max(GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size()) - 1);
273 }
274
275 this->exx_objs.clear();
276 this->coulomb_settings.clear();
277 this->coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2]
278 = std::make_pair(true,
279 ExxLriDetail::build_center2_cut_coulomb_param(
280 this->info.coulomb_param, ucell, kv_in));
281
282 this->MGT = std::make_shared<ORB_gaunt_table>();
283 const auto center2_settings = this->coulomb_settings.find(Conv_Coulomb_Pot_K::Coulomb_Method::Center2);
284 if (center2_settings == this->coulomb_settings.end())
285 {
286 throw std::invalid_argument("Exx_LRI::init_spencer failed to prepare Center2 settings.");
287 }
288 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp_spencer(
289 this->abfs,
290 center2_settings->second.second,
291 this->info.ccp_rmesh_times);
292 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.set_orbitals(
293 ucell,
294 orb,
295 this->lcaos,
296 this->abfs,
297 this->exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].abfs_ccp,
298 this->info.kmesh_times,
299 this->MGT,
300 center2_settings->second.first);
301
302 ModuleBase::timer::end("Exx_LRI", "init_spencer");
303}
304
305template<typename Tdata>
306void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
307 const bool write_cv)
308{
309 ModuleBase::TITLE("Exx_LRI","cal_exx_ions");
310 ModuleBase::timer::start("Exx_LRI", "cal_exx_ions");
311
312 // init_radial_table_ions( cal_atom_centres_core(atom_pairs_core_origin), atom_pairs_core_origin );
313
314 // this->m_abfsabfs.init_radial_table(Rradial);
315 // this->m_abfslcaos_lcaos.init_radial_table(Rradial);
316
317 std::vector<TA> atoms(ucell.nat);
318 for(int iat=0; iat<ucell.nat; ++iat)
319 { atoms[iat] = iat; }
320 std::map<TA,TatomR> atoms_pos;
321 for(int iat=0; iat<ucell.nat; ++iat)
322 { atoms_pos[iat] = RI_Util::Vector3_to_array3( ucell.atoms[ ucell.iat2it[iat] ].tau[ ucell.iat2ia[iat] ] ); }
323 const std::array<TatomR,Ndim> latvec
324 = {RI_Util::Vector3_to_array3(ucell.a1),
325 RI_Util::Vector3_to_array3(ucell.a2),
326 RI_Util::Vector3_to_array3(ucell.a3)};
327 const std::array<Tcell,Ndim> period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]};
328
329 this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
330
331 // std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour.
332 const std::array<Tcell,Ndim> period_Vs = LRI_CV_Tools::cal_latvec_range<Tcell>(1+this->info.ccp_rmesh_times, ucell, orb_cutoff_);
333 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
334 list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
335
336 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Vs;
337 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>> dVs;
338 for(const auto &settings_list : this->coulomb_settings)
339 {
340 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
341 Vs_temp = this->exx_objs[settings_list.first].cv.cal_Vs(ucell,
342 list_As_Vs.first, list_As_Vs.second[0],
343 {{"writable_Vws",true}});
344 this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs_temp);
345 if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald)
346 {
347 this->exx_objs[settings_list.first].evq.init_ions(ucell, period_Vs);
348 const auto &coulomb_param = settings_list.second.second;
349 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald;
350 for(const auto &param_list : coulomb_param)
351 {
352 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald_temp;
353 switch(param_list.first)
354 {
355 case Conv_Coulomb_Pot_K::Coulomb_Type::Fock:
356 {
357 double chi = this->exx_objs[settings_list.first].evq.get_singular_chi(ucell, param_list.second, 2.0);
358 Vs_ewald_temp = this->exx_objs[settings_list.first].evq.cal_Vs(ucell, chi, Vs_temp);
359 break;
360 }
361 default:
362 {
363 throw std::invalid_argument( std::string(__FILE__) + " line " + std::to_string(__LINE__) );
364 }
365 }
366 // Vs_temp cannot be covered here
367 Vs_ewald = Vs_ewald.empty() ? Vs_ewald_temp : LRI_CV_Tools::add(Vs_ewald, Vs_ewald_temp);
368 }
369 Vs_temp = Vs_ewald;
370 }
371 Vs = Vs.empty() ? Vs_temp : LRI_CV_Tools::add(Vs, Vs_temp);
372
373 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
374 {
375 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>>
376 dVs_temp = this->exx_objs[settings_list.first].cv.cal_dVs(ucell,
377 list_As_Vs.first, list_As_Vs.second[0],
378 {{"writable_dVws",true}});
379 this->exx_objs[settings_list.first].cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs_temp);
380 dVs = dVs.empty() ? dVs_temp : LRI_CV_Tools::add(dVs, dVs_temp);
381 }
382 }
383 if (write_cv && GlobalV::MY_RANK == 0)
384 { LRI_CV_Tools::write_Vs_abf(Vs, PARAM.globalv.global_out_dir + "Vs"); }
385 this->exx_lri.set_Vs(std::move(Vs), this->info.V_threshold);
386
387 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
388 {
389 std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
390 dVs_order = LRI_CV_Tools::change_order(std::move(dVs));
391 this->exx_lri.set_dVs(std::move(dVs_order), this->info.V_grad_threshold);
392 if(PARAM.inp.cal_stress)
393 {
394 std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs_order);
395 this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold);
396 }
397 }
398
399 const std::array<Tcell,Ndim> period_Cs = LRI_CV_Tools::cal_latvec_range<Tcell>(2, ucell,orb_cutoff_);
400 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
401 list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
402
403 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Cs;
404 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;
405 for(const auto &settings_list : this->coulomb_settings)
406 {
407 if(settings_list.second.first)
408 {
409 std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
410 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
411 Cs_dCs = this->exx_objs[settings_list.first].cv.cal_Cs_dCs(
412 ucell,
413 list_As_Cs.first, list_As_Cs.second[0],
414 {{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
415 {"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}});
416 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs_temp = std::get<0>(Cs_dCs);
417 this->exx_objs[settings_list.first].cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs_temp);
418 Cs = Cs.empty() ? Cs_temp : LRI_CV_Tools::add(Cs, Cs_temp);
419
420 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
421 {
422 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> &dCs_temp = std::get<1>(Cs_dCs);
423 this->exx_objs[settings_list.first].cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs_temp);
424 dCs = dCs.empty() ? dCs_temp : LRI_CV_Tools::add(dCs, dCs_temp);
425 }
426 }
427 }
428 if (write_cv && GlobalV::MY_RANK == 0)
429 { LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs"); }
430 this->exx_lri.set_Cs(std::move(Cs), this->info.C_threshold);
431
432 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
433 {
434 std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
435 dCs_order = LRI_CV_Tools::change_order(std::move(dCs));
436 this->exx_lri.set_dCs(std::move(dCs_order), this->info.C_grad_threshold);
437 if(PARAM.inp.cal_stress)
438 {
439 std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs_order);
440 this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold);
441 }
442 }
443 ModuleBase::timer::end("Exx_LRI", "cal_exx_ions");
444}
445
446 #if 0
447 template <typename Tdata>
448 void Exx_LRI<Tdata>::cal_cut_coulomb_cs(std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_cut,
449 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs,
450 const UnitCell& ucell,
451 const bool write_cv)
452{
453 ModuleBase::TITLE("Exx_LRI", "cal_cut_coulomb_cs");
454 ModuleBase::timer::start("Exx_LRI", "cal_cut_coulomb_cs");
455
456 std::vector<TA> atoms(ucell.nat);
457 for(int iat=0; iat<ucell.nat; ++iat)
458 atoms[iat] = iat;
459 std::map<TA,TatomR> atoms_pos;
460 for(int iat=0; iat<ucell.nat; ++iat)
461 atoms_pos[iat] = RI_Util::Vector3_to_array3( ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]] );
462 const std::array<TatomR,Ndim> latvec
463 = {RI_Util::Vector3_to_array3(ucell.a1),
464 RI_Util::Vector3_to_array3(ucell.a2),
465 RI_Util::Vector3_to_array3(ucell.a3)};
466 const std::array<Tcell,Ndim> period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]};
467
468 this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
469
470 // std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour.
471 const std::array<Tcell, Ndim> period_Vs
472 = LRI_CV_Tools::cal_latvec_range<Tcell>(1 + this->info.ccp_rmesh_times, ucell, orb_cutoff_);
473 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
474 list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
475
476 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>> dVs;
477 for(const auto &settings_list : this->coulomb_settings)
478 {
479 if(!settings_list.second.first) continue;
480 std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
481 Vs_temp = this->exx_objs[settings_list.first].cv.cal_Vs(ucell,
482 list_As_Vs.first, list_As_Vs.second[0],
483 {{"writable_Vws",true}});
484 this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs_temp);
485
486 // ======rotate ABFs begin======
487 int flag = 0;
488 for (const auto& IJRc: this->exx_objs[settings_list.first].cv.Vws)
489 {
490 const TA& I = IJRc.first;
491 const auto& JRc = IJRc.second;
492 for (const auto& JRc_tensor: JRc)
493 {
494 const TA& J = JRc_tensor.first;
495 const auto Rc = JRc_tensor.second;
496 for (const auto& Rc_tensor: Rc)
497 {
498 const auto& R = Rc_tensor.first;
499 flag += 1;
500 }
501 }
502 }
503 std::cout << "Coulomb: number of atom-pairs inside atomic overlap is " << flag << ". " << std::endl;
504 if (this->info.coul_moment == true)
505 {
506 double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0);
507 // To cal Cs, we still cal all Vs(R) in r space
508 // moment_abfs->cal_VR(ucell,
509 // this->abfs,
510 // list_As_Vs,
511 // orb_cutoff_,
512 // hf_Rcut,
513 // this->exx_objs[settings_list.first].cv,
514 // Vs_cut);
515 delete moment_abfs;
516 moment_abfs = nullptr;
517 malloc_trim(0);
518 }
519
520 flag = 0;
521 for (const auto& IJRc: this->exx_objs[settings_list.first].cv.Vws)
522 {
523 const auto& JRc = IJRc.second;
524 for (const auto& JRc_tensor: JRc)
525 {
526 const auto Rc = JRc_tensor.second;
527 for (const auto& Rc_tensor: Rc)
528 {
529 flag += 1;
530 }
531 }
532 }
533 std::cout << "Coulomb: number of all atom-pairs is " << flag << ". " << std::endl;
534 // ======rotate ABFs end======
535
536 Vs_cut = Vs.empty() ? Vs_temp : LRI_CV_Tools::add(Vs_cut, Vs_temp);
537
538 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
539 {
540 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>>
541 dVs_temp = this->exx_objs[settings_list.first].cv.cal_dVs(ucell,
542 list_As_Vs.first, list_As_Vs.second[0],
543 {{"writable_dVws",true}});
544 this->exx_objs[settings_list.first].cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs_temp);
545 dVs = dVs.empty() ? dVs_temp : LRI_CV_Tools::add(dVs, dVs_temp);
546 }
547 }
548
549 if (write_cv && GlobalV::MY_RANK == 0)
550 {
551 LRI_CV_Tools::write_Vs_abf(Vs_cut, PARAM.globalv.global_out_dir + "Vs_cut");
552 }
553 this->exx_lri.set_Vs(std::move(Vs_cut), this->info.V_threshold);
554
555 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
556 {
557 std::map<TA,std::map<TAC,std::array<RI::Tensor<Tdata>,Ndim>>> dVs
558 = this->exx_objs[coulomb_method].cv.cal_dVs(ucell,
559 list_As_Vs.first, list_As_Vs.second[0],
560 {{"writable_dVws",true}});
561 this->exx_objs[coulomb_method].cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs);
562
563 std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,Ndim> dVs_order
564 = LRI_CV_Tools::change_order(std::move(dVs));
565 this->exx_lri.set_dVs(std::move(dVs_order), this->info.V_grad_threshold);
566 if(PARAM.inp.cal_stress)
567 {
568 std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs
569 = LRI_CV_Tools::cal_dMRs(ucell,dVs_order);
570 this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold);
571 }
572 }
573
574 const std::array<Tcell,Ndim> period_Cs = LRI_CV_Tools::cal_latvec_range<Tcell>(2, ucell,orb_cutoff_);
575 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
576 list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
577 std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,
578 std::map<TA,std::map<TAC,std::array<RI::Tensor<Tdata>,3>>>>
579 Cs_dCs = this->exx_objs[coulomb_method].cv.cal_Cs_dCs(ucell,
580 list_As_Cs.first, list_As_Cs.second[0],
581 {{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
582 {"writable_Cws",true},
583 {"writable_dCws",true},
584 {"writable_Vws",false},
585 {"writable_dVws",false}});
586 Cs = std::get<0>(Cs_dCs);
587 this->exx_objs[coulomb_method].cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs);
588 if(write_cv && GlobalV::MY_RANK==0)
589 {
590 LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs");
591 }
592 this->exx_lri.set_Cs(Cs, this->info.C_threshold);
593
594 if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
595 {
596 std::map<TA,std::map<TAC,std::array<RI::Tensor<Tdata>,3>>>& dCs = std::get<1>(Cs_dCs);
597 this->exx_objs[coulomb_method].cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs);
598 std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,Ndim> dCs_order
599 = LRI_CV_Tools::change_order(std::move(dCs));
600 this->exx_lri.set_dCs(std::move(dCs_order), this->info.C_grad_threshold);
601 if(PARAM.inp.cal_stress)
602 {
603 std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs
604 = LRI_CV_Tools::cal_dMRs(ucell,dCs_order);
605 this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold);
606 }
607 }
608 ModuleBase::timer::end("Exx_LRI", "cal_cut_coulomb_cs");
609}
610
611template <typename Tdata>
612void Exx_LRI<Tdata>::cal_ewald_coulomb(std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_full,
613 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs,
614 const UnitCell& ucell,
615 const bool write_cv)
616{
617 ModuleBase::TITLE("Exx_LRI", "cal_ewald_coulomb");
618 ModuleBase::timer::start("Exx_LRI", "cal_ewald_coulomb");
619
620 std::vector<TA> atoms(ucell.nat);
621 for (int iat = 0; iat < ucell.nat; ++iat)
622 atoms[iat] = iat;
623 std::map<TA, TatomR> atoms_pos;
624 for (int iat = 0; iat < ucell.nat; ++iat)
625 atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
626 const std::array<TatomR, Ndim> latvec = {RI_Util::Vector3_to_array3(ucell.a1),
627 RI_Util::Vector3_to_array3(ucell.a2),
628 RI_Util::Vector3_to_array3(ucell.a3)};
629 const std::array<Tcell, Ndim> period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]};
630
631 this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
632
633 // std::max(3) for gamma_only, list_A2 should contain cell {-1,0,1}. In the future distribute will be neighbour.
634 const std::array<Tcell, Ndim> period_Vs
635 = LRI_CV_Tools::cal_latvec_range<Tcell>(1 + this->info.ccp_rmesh_times, ucell, orb_cutoff_);
636 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA, std::array<Tcell, Ndim>>>>> list_As_Vs
637 = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
638
639 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>> dVs;
640 for (const auto& settings_list: this->coulomb_settings)
641 {
642 if (!settings_list.second.first)
643 continue;
644 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_temp
645 = this->exx_objs[settings_list.first].cv.cal_Vs(ucell,
646 list_As_Vs.first,
647 list_As_Vs.second[0],
648 {{"writable_Vws", true}});
649 this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell, Vs_temp);
650
651 // ======rotate ABFs begin======
652 int flag = 0;
653 for (const auto& IJRc: this->exx_objs[settings_list.first].cv.Vws)
654 {
655 const TA& I = IJRc.first;
656 const auto& JRc = IJRc.second;
657 for (const auto& JRc_tensor: JRc)
658 {
659 const TA& J = JRc_tensor.first;
660 const auto Rc = JRc_tensor.second;
661 for (const auto& Rc_tensor: Rc)
662 {
663 const auto& R = Rc_tensor.first;
664 flag += 1;
665 }
666 }
667 }
668 std::cout << "Coulomb: number of atom-pairs inside atomic overlap is " << flag << ". " << std::endl;
669 if (this->info.coul_moment == true)
670 {
671 double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0);
672 // To cal Cs, we still cal all Vs(R) in r space
673 // moment_abfs->cal_VR(ucell,
674 // this->abfs,
675 // list_As_Vs,
676 // orb_cutoff_,
677 // hf_Rcut,
678 // this->exx_objs[settings_list.first].cv,
679 // Vs_full);
680 delete moment_abfs;
681 moment_abfs = nullptr;
682 malloc_trim(0);
683 }
684
685 flag = 0;
686 for (const auto& IJRc: this->exx_objs[settings_list.first].cv.Vws)
687 {
688 const auto& JRc = IJRc.second;
689 for (const auto& JRc_tensor: JRc)
690 {
691 const auto Rc = JRc_tensor.second;
692 for (const auto& Rc_tensor: Rc)
693 {
694 flag += 1;
695 }
696 }
697 }
698 std::cout << "Coulomb: number of all atom-pairs is " << flag << ". " << std::endl;
699 // ======rotate ABFs end======
700
701 if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald)
702 {
703 this->exx_objs[settings_list.first].evq.init_ions(ucell, period_Vs);
704 const auto& coulomb_param = settings_list.second.second;
705 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald;
706 for (const auto& param_list: coulomb_param)
707 {
708 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald_temp;
709 switch (param_list.first)
710 {
711 case Conv_Coulomb_Pot_K::Coulomb_Type::Fock: {
712 double chi
713 = this->exx_objs[settings_list.first].evq.get_singular_chi(ucell, param_list.second, 2.0);
714 Vs_ewald_temp = this->exx_objs[settings_list.first].evq.cal_Vs(ucell, chi, Vs_temp);
715 break;
716 }
717 default: {
718 throw std::invalid_argument(std::string(__FILE__) + " line " + std::to_string(__LINE__));
719 }
720 }
721 // Vs_temp cannot be covered here
722 Vs_ewald = Vs_ewald.empty() ? Vs_ewald_temp : LRI_CV_Tools::add(Vs_ewald, Vs_ewald_temp);
723 }
724 Vs_temp = Vs_ewald;
725 }
726
727 Vs_full = Vs.empty() ? Vs_temp : LRI_CV_Tools::add(Vs_full, Vs_temp);
728 }
729
730 if (write_cv && GlobalV::MY_RANK == 0)
731 {
732 LRI_CV_Tools::write_Vs_abf(Vs_full, PARAM.globalv.global_out_dir + "Vs_full");
733 }
734 // this->exx_lri.set_Vs(std::move(Vs_full), this->info.V_threshold);
735
736 // const std::array<Tcell,Ndim> period_Cs = LRI_CV_Tools::cal_latvec_range<Tcell>(2, ucell,orb_cutoff_);
737 // const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
738 // list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
739 // std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,
740 // std::map<TA,std::map<TAC,std::array<RI::Tensor<Tdata>,3>>>>
741 // Cs_dCs = this->exx_objs[coulomb_method].cv.cal_Cs_dCs(ucell,
742 // list_As_Cs.first, list_As_Cs.second[0],
743 // {{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
744 // {"writable_Cws",true},
745 // {"writable_dCws",true},
746 // {"writable_Vws",false},
747 // {"writable_dVws",false}});
748 // Cs = std::get<0>(Cs_dCs);
749 // this->exx_objs[coulomb_method].cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs);
750 // if(write_cv && GlobalV::MY_RANK==0)
751 // {
752 // LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs");
753 // }
754 // this->exx_lri.set_Cs(Cs, this->info.C_threshold);
755 ModuleBase::timer::end("Exx_LRI", "cal_ewald_coulomb");
756}
757 #endif
758
759template <typename Tdata>
760void Exx_LRI<Tdata>::cal_cut_coulomb_cs(std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_cut,
761 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs,
762 const UnitCell& ucell,
763 const bool write_cv)
764{
765 ModuleBase::TITLE("Exx_LRI", "cal_cut_coulomb_cs");
766 ModuleBase::timer::start("Exx_LRI", "cal_cut_coulomb_cs");
767
768 std::vector<TA> atoms(ucell.nat);
769 for (int iat = 0; iat < ucell.nat; ++iat)
770 {
771 atoms[iat] = iat;
772 }
773 std::map<TA, TatomR> atoms_pos;
774 for (int iat = 0; iat < ucell.nat; ++iat)
775 {
776 atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
777 }
778 const std::array<TatomR, Ndim> latvec = {RI_Util::Vector3_to_array3(ucell.a1),
779 RI_Util::Vector3_to_array3(ucell.a2),
780 RI_Util::Vector3_to_array3(ucell.a3)};
781 const std::array<Tcell, Ndim> period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]};
782
783 this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
784
785 const auto center2_method = Conv_Coulomb_Pot_K::Coulomb_Method::Center2;
786 auto center2_obj_it = this->exx_objs.find(center2_method);
787 if (center2_obj_it == this->exx_objs.end())
788 {
789 throw std::invalid_argument(std::string(__FILE__) + " line " + std::to_string(__LINE__));
790 }
791
792 const std::array<Tcell, Ndim> period_Vs
793 = LRI_CV_Tools::cal_latvec_range<Tcell>(1 + this->info.ccp_rmesh_times, ucell, orb_cutoff_);
794 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA, std::array<Tcell, Ndim>>>>> list_As_Vs
795 = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
796
797 Vs_cut = center2_obj_it->second.cv.cal_Vs(
798 ucell,
799 list_As_Vs.first,
800 list_As_Vs.second[0],
801 {{"writable_Vws", true}});
802 center2_obj_it->second.cv.Vws = LRI_CV_Tools::get_CVws(ucell, Vs_cut);
803 if (write_cv && GlobalV::MY_RANK == 0)
804 {
805 LRI_CV_Tools::write_Vs_abf(Vs_cut, PARAM.globalv.global_out_dir + "Vs_cut");
806 }
807 this->exx_lri.set_Vs(Vs_cut, this->info.V_threshold);
808
809 const std::array<Tcell, Ndim> period_Cs = LRI_CV_Tools::cal_latvec_range<Tcell>(2, ucell, orb_cutoff_);
810 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA, std::array<Tcell, Ndim>>>>> list_As_Cs
811 = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
812 std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
813 std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
814 Cs_dCs = center2_obj_it->second.cv.cal_Cs_dCs(
815 ucell,
816 list_As_Cs.first,
817 list_As_Cs.second[0],
818 {{"cal_dC", false},
819 {"writable_Cws", true},
820 {"writable_dCws", true},
821 {"writable_Vws", false},
822 {"writable_dVws", false}});
823 Cs = std::get<0>(Cs_dCs);
824 center2_obj_it->second.cv.Cws = LRI_CV_Tools::get_CVws(ucell, Cs);
825 if (write_cv && GlobalV::MY_RANK == 0)
826 {
827 LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs");
828 }
829 this->exx_lri.set_Cs(Cs, this->info.C_threshold);
830
831 ModuleBase::timer::end("Exx_LRI", "cal_cut_coulomb_cs");
832}
833
834template <typename Tdata>
835void Exx_LRI<Tdata>::cal_ewald_coulomb(std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Vs_full,
836 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>& Cs,
837 const UnitCell& ucell,
838 const bool write_cv)
839{
840 ModuleBase::TITLE("Exx_LRI", "cal_ewald_coulomb");
841 ModuleBase::timer::start("Exx_LRI", "cal_ewald_coulomb");
842
843 std::vector<TA> atoms(ucell.nat);
844 for (int iat = 0; iat < ucell.nat; ++iat)
845 {
846 atoms[iat] = iat;
847 }
848 std::map<TA, TatomR> atoms_pos;
849 for (int iat = 0; iat < ucell.nat; ++iat)
850 {
851 atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
852 }
853 const std::array<TatomR, Ndim> latvec = {RI_Util::Vector3_to_array3(ucell.a1),
854 RI_Util::Vector3_to_array3(ucell.a2),
855 RI_Util::Vector3_to_array3(ucell.a3)};
856 const std::array<Tcell, Ndim> period = {this->p_kv->nmp[0], this->p_kv->nmp[1], this->p_kv->nmp[2]};
857
858 this->exx_lri.set_parallel(this->mpi_comm, atoms_pos, latvec, period);
859
860 const std::array<Tcell, Ndim> period_Vs
861 = LRI_CV_Tools::cal_latvec_range<Tcell>(1 + this->info.ccp_rmesh_times, ucell, orb_cutoff_);
862 const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA, std::array<Tcell, Ndim>>>>> list_As_Vs
863 = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
864
865 for (const auto& settings_list : this->coulomb_settings)
866 {
867 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_temp
868 = this->exx_objs[settings_list.first].cv.cal_Vs(
869 ucell,
870 list_As_Vs.first,
871 list_As_Vs.second[0],
872 {{"writable_Vws", true}});
873 this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell, Vs_temp);
874
875 if (settings_list.first == Conv_Coulomb_Pot_K::Coulomb_Method::Ewald)
876 {
877 this->exx_objs[settings_list.first].evq.init_ions(ucell, period_Vs);
878 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald;
879 for (const auto& param_list : settings_list.second.second)
880 {
881 std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs_ewald_temp;
882 switch (param_list.first)
883 {
884 case Conv_Coulomb_Pot_K::Coulomb_Type::Fock:
885 {
886 double chi = this->exx_objs[settings_list.first].evq.get_singular_chi(ucell, param_list.second, 2.0);
887 Vs_ewald_temp = this->exx_objs[settings_list.first].evq.cal_Vs(ucell, chi, Vs_temp);
888 break;
889 }
890 default:
891 {
892 throw std::invalid_argument(std::string(__FILE__) + " line " + std::to_string(__LINE__));
893 }
894 }
895 Vs_ewald = Vs_ewald.empty() ? Vs_ewald_temp : LRI_CV_Tools::add(Vs_ewald, Vs_ewald_temp);
896 }
897 Vs_temp = Vs_ewald;
898 }
899
900 Vs_full = Vs_full.empty() ? Vs_temp : LRI_CV_Tools::add(Vs_full, Vs_temp);
901 }
902
903 if (write_cv && GlobalV::MY_RANK == 0)
904 {
905 LRI_CV_Tools::write_Vs_abf(Vs_full, PARAM.globalv.global_out_dir + "Vs_full");
906 }
907
908 Cs.clear();
909 ModuleBase::timer::end("Exx_LRI", "cal_ewald_coulomb");
910}
911
912template<typename Tdata>
913void Exx_LRI<Tdata>::cal_exx_elec(const std::vector<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>>& Ds,
914 const UnitCell& ucell,
915 const Parallel_Orbitals& pv,
916 const ModuleSymmetry::Symmetry_rotation* p_symrot)
917{
918 ModuleBase::TITLE("Exx_LRI","cal_exx_elec");
919 ModuleBase::timer::start("Exx_LRI", "cal_exx_elec");
920
921 const std::vector<std::tuple<std::set<TA>, std::set<TA>>> judge = RI_2D_Comm::get_2D_judge(ucell,pv);
922
923 if(p_symrot)
924 { this->exx_lri.set_symmetry(true, p_symrot->get_irreducible_sector()); }
925 else
926 { this->exx_lri.set_symmetry(false, {}); }
927
928 this->Hexxs.resize(PARAM.inp.nspin);
929 this->Eexx = 0;
930 for(int is=0; is<PARAM.inp.nspin; ++is)
931 {
932 const std::string suffix = ((PARAM.inp.cal_force || PARAM.inp.cal_stress) ? std::to_string(is) : "");
933
934 this->exx_lri.set_Ds(Ds[is], this->info.dm_threshold, suffix);
935 this->exx_lri.cal_Hs({ "","",suffix });
936
937 if (!p_symrot)
938 {
939 this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first(
940 this->mpi_comm, std::move(this->exx_lri.Hs), std::get<0>(judge[is]), std::get<1>(judge[is]));
941 }
942 else
943 {
944 // reduce but not repeat
945 auto Hs_a2D = this->exx_lri.post_2D.set_tensors_map2(this->exx_lri.Hs);
946 // rotate locally without repeat
947 Hs_a2D = p_symrot->restore_HR(ucell.symm, ucell.atoms, ucell.st, 'H', Hs_a2D);
948 // cal energy using full Hs without repeat
949 this->exx_lri.energy = this->exx_lri.post_2D.cal_energy(
950 this->exx_lri.post_2D.saves["Ds_" + suffix],
951 this->exx_lri.post_2D.set_tensors_map2(Hs_a2D));
952 // get repeated full Hs for abacus
953 this->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first(
954 this->mpi_comm, std::move(Hs_a2D), std::get<0>(judge[is]), std::get<1>(judge[is]));
955 }
956 this->Eexx += std::real(this->exx_lri.energy);
957 post_process_Hexx(this->Hexxs[is]);
958 }
959 this->Eexx = post_process_Eexx(this->Eexx);
960 this->exx_lri.set_symmetry(false, {});
961 ModuleBase::timer::end("Exx_LRI", "cal_exx_elec");
962}
963
964template<typename Tdata>
965void Exx_LRI<Tdata>::post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> &Hexxs_io ) const
966{
967 ModuleBase::TITLE("Exx_LRI","post_process_Hexx");
968 constexpr Tdata frac = -1 * 2; // why? Hartree to Ry?
969 const std::function<void(RI::Tensor<Tdata>&)>
970 multiply_frac = [&frac](RI::Tensor<Tdata> &t)
971 { t = t*frac; };
972 RI::Map_Operator::for_each( Hexxs_io, multiply_frac );
973}
974
975template<typename Tdata>
976double Exx_LRI<Tdata>::post_process_Eexx(const double& Eexx_in) const
977{
978 ModuleBase::TITLE("Exx_LRI","post_process_Eexx");
979 const double SPIN_multiple = std::map<int, double>{ {1,2}, {2,1}, {4,1} }.at(PARAM.inp.nspin); // why?
980 const double frac = -SPIN_multiple;
981 return frac * Eexx_in;
982}
983
984/*
985post_process_old
986{
987 // D
988 const std::map<int,double> SPIN_multiple = {{1,0.5}, {2,1}, {4,1}}; // ???
989 DR *= SPIN_multiple.at(NSPIN);
990
991 // H
992 HR *= -2;
993
994 // E
995 const std::map<int,double> SPIN_multiple = {{1,2}, {2,1}, {4,1}}; // ???
996 energy *= SPIN_multiple.at(PARAM.inp.nspin); // ?
997 energy /= 2; // /2 for Ry
998}
999*/
1000
1001template<typename Tdata>
1002void Exx_LRI<Tdata>::cal_exx_force(const int& nat)
1003{
1004 ModuleBase::TITLE("Exx_LRI","cal_exx_force");
1005 ModuleBase::timer::start("Exx_LRI", "cal_exx_force");
1006
1007 this->force_exx.create(nat, Ndim);
1008 for(int is=0; is<PARAM.inp.nspin; ++is)
1009 {
1010 this->exx_lri.cal_force({"","",std::to_string(is),"",""});
1011 for(std::size_t idim=0; idim<Ndim; ++idim) {
1012 for(const auto &force_item : this->exx_lri.force[idim]) {
1013 this->force_exx(force_item.first, idim) += std::real(force_item.second);
1014 } }
1015 }
1016
1017 const double SPIN_multiple = std::map<int,double>{{1,2}, {2,1}, {4,1}}.at(PARAM.inp.nspin); // why?
1018 const double frac = -2 * SPIN_multiple; // why?
1019 this->force_exx *= frac;
1020 ModuleBase::timer::end("Exx_LRI", "cal_exx_force");
1021}
1022
1023
1024template<typename Tdata>
1025void Exx_LRI<Tdata>::cal_exx_stress(const double& omega, const double& lat0)
1026{
1027 ModuleBase::TITLE("Exx_LRI","cal_exx_stress");
1028 ModuleBase::timer::start("Exx_LRI", "cal_exx_stress");
1029
1030 this->stress_exx.create(Ndim, Ndim);
1031 for(int is=0; is<PARAM.inp.nspin; ++is)
1032 {
1033 this->exx_lri.cal_stress({"","",std::to_string(is),"",""});
1034 for(std::size_t idim0=0; idim0<Ndim; ++idim0) {
1035 for(std::size_t idim1=0; idim1<Ndim; ++idim1) {
1036 this->stress_exx(idim0,idim1) += std::real(this->exx_lri.stress(idim0,idim1));
1037 } }
1038 }
1039
1040 const double SPIN_multiple = std::map<int,double>{{1,2}, {2,1}, {4,1}}.at(PARAM.inp.nspin); // why?
1041 const double frac = 2 * SPIN_multiple / omega * lat0; // why?
1042 this->stress_exx *= frac;
1043
1044 ModuleBase::timer::end("Exx_LRI", "cal_exx_stress");
1045}
1046
1047/*
1048template<typename Tdata>
1049std::vector<std::vector<int>> Exx_LRI<Tdata>::get_abfs_nchis() const
1050{
1051 std::vector<std::vector<int>> abfs_nchis;
1052 for (const auto& abfs_T : this->abfs)
1053 {
1054 std::vector<int> abfs_nchi_T;
1055 for (const auto& abfs_L : abfs_T)
1056 { abfs_nchi_T.push_back(abfs_L.size()); }
1057 abfs_nchis.push_back(abfs_nchi_T);
1058 }
1059 return abfs_nchis;
1060}
1061*/
1062
1063#endif
Exx_LRI.h
Mix_DMk_2D.h
RI_2D_Comm.h
RI_Util.h
Atom::tau
std::vector< ModuleBase::Vector3< double > > tau
Definition atom_spec.h:35
Exx_Abfs::Construct_Orbs::print_orbs_size
static void print_orbs_size(const UnitCell &ucell, const std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > &orbs, std::ostream &os)
Definition exx_abfs-construct_orbs.cpp:475
Exx_Abfs::Construct_Orbs::filter_empty_orbs
static void filter_empty_orbs(std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > &orbs)
Definition exx_abfs-construct_orbs.cpp:545
Exx_Abfs::Construct_Orbs::abfs_same_atom
static std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > abfs_same_atom(const UnitCell &ucell, const LCAO_Orbitals &orb, const std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > &lcaos, const double kmesh_times_mot, const double times_threshold=0)
Definition exx_abfs-construct_orbs.cpp:80
Exx_Abfs::Construct_Orbs::change_orbs
static std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > change_orbs(const LCAO_Orbitals &orb_in, const double kmesh_times)
Definition exx_abfs-construct_orbs.cpp:10
Exx_Abfs::IO::construct_abfs
static std::vector< std::vector< std::vector< Numerical_Orbital_Lm > > > construct_abfs(const LCAO_Orbitals &orbs, const std::vector< std::string > &files_abfs, const double kmesh_times=1)
Definition exx_abfs-io.cpp:12
Exx_LRI::post_process_Hexx
void post_process_Hexx(std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > &Hexxs_io) const
Definition Exx_LRI.hpp:965
Exx_LRI::init
void init(const MPI_Comm &mpi_comm_in, const UnitCell &ucell, const K_Vectors &kv_in, const LCAO_Orbitals &orb)
Definition Exx_LRI.hpp:85
Exx_LRI::cal_cut_coulomb_cs
void cal_cut_coulomb_cs(std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > &Vs_cut_IJR, std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > &Cs, const UnitCell &ucell, const bool write_cv=false)
Definition Exx_LRI.hpp:760
Exx_LRI::init_spencer
void init_spencer(const MPI_Comm &mpi_comm_in, const UnitCell &ucell, const K_Vectors &kv_in, const LCAO_Orbitals &orb)
Definition Exx_LRI.hpp:179
Exx_LRI::cal_exx_elec
void cal_exx_elec(const std::vector< std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > > &Ds, const UnitCell &ucell, const Parallel_Orbitals &pv, const ModuleSymmetry::Symmetry_rotation *p_symrot=nullptr)
Definition Exx_LRI.hpp:913
Exx_LRI::post_process_Eexx
double post_process_Eexx(const double &Eexx_in) const
Definition Exx_LRI.hpp:976
Exx_LRI::TAC
std::pair< TA, TC > TAC
Definition Exx_LRI.h:59
Exx_LRI::TA
int TA
Definition Exx_LRI.h:55
Exx_LRI::cal_exx_ions
void cal_exx_ions(const UnitCell &ucell, const bool write_cv=false)
Definition Exx_LRI.hpp:306
Exx_LRI::cal_exx_force
void cal_exx_force(const int &nat)
Definition Exx_LRI.hpp:1002
Exx_LRI::cal_ewald_coulomb
void cal_ewald_coulomb(std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > &Vs_full_IJR, std::map< TA, std::map< TAC, RI::Tensor< Tdata > > > &Cs, const UnitCell &ucell, const bool write_cv=false)
Definition Exx_LRI.hpp:835
Exx_LRI::cal_exx_stress
void cal_exx_stress(const double &omega, const double &lat0)
Definition Exx_LRI.hpp:1025
K_Vectors
Definition klist.h:12
K_Vectors::get_nkstot_full
int get_nkstot_full() const
Definition klist.h:78
LCAO_Orbitals
Definition ORB_read.h:18
LCAO_Orbitals::cutoffs
std::vector< double > cutoffs() const
Definition ORB_read.cpp:39
ModuleBase::timer::start
static void start(void)
Start total time calculation.
Definition timer.cpp:44
ModuleBase::timer::end
static void end(const std::string &class_name_in, const std::string &name_in)
Definition timer.cpp:109
ModuleSymmetry::Symmetry_rotation
Definition symmetry_rotation.h:16
ModuleSymmetry::Symmetry_rotation::restore_HR
std::map< int, std::map< std::pair< int, TC >, RI::Tensor< Tdata > > > restore_HR(const Symmetry &symm, const Atom *atoms, const Statistics &st, const char mode, const std::map< int, std::map< std::pair< int, TC >, RI::Tensor< Tdata > > > &HR_irreduceble) const
Definition symmetry_rotation_R.hpp:40
ModuleSymmetry::Symmetry_rotation::get_irreducible_sector
const std::map< Tap, std::set< TC > > & get_irreducible_sector() const
Definition symmetry_rotation.h:23
Parallel_Orbitals
Definition parallel_orbitals.h:9
Parameter::inp
const Input_para & inp
Definition parameter.h:26
Parameter::globalv
const System_para & globalv
Definition parameter.h:30
UnitCell
Definition unitcell.h:15
UnitCell::iat2it
int *& iat2it
Definition unitcell.h:75
UnitCell::atoms
Atom * atoms
Definition unitcell.h:45
UnitCell::symm
ModuleSymmetry::Symmetry symm
Definition unitcell.h:83
UnitCell::nat
int & nat
Definition unitcell.h:74
UnitCell::omega
double & omega
Definition unitcell.h:60
UnitCell::a2
ModuleBase::Vector3< double > & a2
Definition unitcell.h:64
UnitCell::a3
ModuleBase::Vector3< double > & a3
Definition unitcell.h:64
UnitCell::iat2ia
int *& iat2ia
Definition unitcell.h:76
UnitCell::a1
ModuleBase::Vector3< double > & a1
Definition unitcell.h:64
UnitCell::st
Statistics st
Definition unitcell.h:72
conv_coulomb_pot_k.h
T
#define T
Definition exp.cpp:237
exx_abfs-construct_orbs.h
exx_abfs-io.h
param
Parameter param
Definition io_system_variable_test.cpp:38
Conv_Coulomb_Pot_K::Coulomb_Type::Fock
@ Fock
Conv_Coulomb_Pot_K::Coulomb_Method::Ewald
@ Ewald
Conv_Coulomb_Pot_K::Coulomb_Method::Center2
@ Center2
Conv_Coulomb_Pot_K::cal_orbs_ccp_spencer
T cal_orbs_ccp_spencer(const T &orbs, const std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > &coulomb_param, const double rmesh_times)
Conv_Coulomb_Pot_K::cal_orbs_ccp
T cal_orbs_ccp(const T &orbs, const std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > &coulomb_param, const double rmesh_times)
ExxLriDetail
Definition Exx_LRI.hpp:34
ExxLriDetail::CoulombParam
std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > CoulombParam
Definition Exx_LRI.hpp:36
ExxLriDetail::trim_malloc_cache
void trim_malloc_cache()
Definition Exx_LRI.hpp:38
ExxLriDetail::default_spencer_rcut
double default_spencer_rcut(const UnitCell &ucell, const K_Vectors &kv)
Definition Exx_LRI.hpp:45
ExxLriDetail::build_center2_cut_coulomb_param
CoulombParam build_center2_cut_coulomb_param(const CoulombParam &coulomb_param, const UnitCell &ucell, const K_Vectors &kv, bool *synthesized_rcut=nullptr)
Definition Exx_LRI.hpp:50
GlobalC::exx_info
Exx_Info exx_info
Definition exx_info.cpp:8
GlobalV::MY_RANK
int MY_RANK
global index of process
Definition global_variable.cpp:21
GlobalV::ofs_running
std::ofstream ofs_running
Definition global_variable.cpp:38
LRI_CV_Tools::add
std::vector< std::array< T, N > > add(const std::vector< std::array< T, N > > &v1, const std::vector< std::array< T, N > > &v2)
Definition LRI_CV_Tools.hpp:222
LRI_CV_Tools::get_CVws
std::map< int, std::map< int, std::map< Abfs::Vector3_Order< double >, RI::Tensor< Tdata > > > > get_CVws(const UnitCell &ucell, const std::map< TA, std::map< std::pair< TA, std::array< Tcell, 3 > >, RI::Tensor< Tdata > > > &CVs)
Definition LRI_CV_Tools.hpp:385
LRI_CV_Tools::change_order
std::array< std::vector< T >, N > change_order(std::vector< std::array< T, N > > &&ds_in)
Definition LRI_CV_Tools.hpp:305
LRI_CV_Tools::get_dCVws
std::map< int, std::map< int, std::map< Abfs::Vector3_Order< double >, std::array< RI::Tensor< Tdata >, 3 > > > > get_dCVws(const UnitCell &ucell, const std::map< TA, std::map< std::pair< TA, std::array< Tcell, 3 > >, std::array< RI::Tensor< Tdata >, 3 > > > &dCVs)
Definition LRI_CV_Tools.hpp:412
LRI_CV_Tools::write_Cs_ao
void write_Cs_ao(const TLRI< T > &Vs, const std::string &file_path)
Definition write_ri_cv.hpp:50
LRI_CV_Tools::cal_dMRs
std::array< std::array< std::map< TA, std::map< std::pair< TA, TC >, RI::Tensor< Tdata > > >, 3 >, 3 > cal_dMRs(const UnitCell &ucell, const std::array< std::map< TA, std::map< std::pair< TA, TC >, RI::Tensor< Tdata > > >, 3 > &dMs)
Definition LRI_CV_Tools.hpp:523
LRI_CV_Tools::write_Vs_abf
void write_Vs_abf(const TLRI< T > &Vs, const std::string &file_path)
Definition write_ri_cv.hpp:105
ModuleBase::GlobalFunc::TO_STRING
std::string TO_STRING(const T &t, const int n=20)
Definition global_function.h:239
ModuleBase::PI
const double PI
Definition constants.h:19
ModuleBase::TITLE
void TITLE(const std::string &class_name, const std::string &function_name, const bool disable)
Definition tool_title.cpp:18
RI_2D_Comm::TA
int TA
Definition RI_2D_Comm.h:25
RI_2D_Comm::get_2D_judge
std::vector< std::tuple< std::set< TA >, std::set< TA > > > get_2D_judge(const UnitCell &ucell, const Parallel_2D &pv)
Definition RI_2D_Comm.cpp:15
RI_Util::Vector3_to_array3
std::array< Tcell, 3 > Vector3_to_array3(const ModuleBase::Vector3< Tcell > &v)
Definition RI_Util.h:32
RI_Util::update_coulomb_param
std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > update_coulomb_param(const std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > &coulomb_param, const UnitCell &ucell, const K_Vectors *p_kv)
Definition RI_Util.hpp:71
RI_Util::update_coulomb_settings
std::map< Conv_Coulomb_Pot_K::Coulomb_Method, std::pair< bool, std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > > > update_coulomb_settings(const std::map< Conv_Coulomb_Pot_K::Coulomb_Type, std::vector< std::map< std::string, std::string > > > &coulomb_param, const UnitCell &ucell, const K_Vectors *p_kv)
Definition RI_Util.hpp:105
PARAM
Parameter PARAM
Definition parameter.cpp:3
TAC
std::pair< int, TC > TAC
Definition ri_cv_io_test.cpp:10
serialization_cereal.h
parallel_orbitals.h
Exx_Info::Exx_Info_RI::abfs_Lmax
int abfs_Lmax
Definition exx_info.h:79
Exx_Info::info_ri
Exx_Info_RI info_ri
Definition exx_info.h:86
Input_para::cal_force
bool cal_force
calculate the force
Definition input_parameter.h:31
Input_para::nspin
int nspin
LDA ; LSDA ; non-linear spin.
Definition input_parameter.h:88
Input_para::cal_stress
bool cal_stress
calculate the stress
Definition input_parameter.h:32
System_para::global_out_dir
std::string global_out_dir
global output directory
Definition system_parameter.h:42
tool_title.h
Generated by doxygen 1.9.8