ABACUS develop
Atomic-orbital Based Ab-initio Computation at UStc
Loading...
Searching...
No Matches
exx_lip.hpp
Go to the documentation of this file.
1//==========================================================
2// AUTHOR : Peize Lin
3
4// DATE : 2015-03-10
5//==========================================================
6
7#ifndef EXX_LIP_HPP
8#define EXX_LIP_HPP
9
10#include "exx_lip.h"
11#include "source_base/vector3.h"
12#include "source_base/global_function.h"
13#include "source_base/vector3.h"
14#include "source_pw/module_pwdft/global.h"
15#include "source_cell/klist.h"
16#include "source_lcao/wavefunc_in_pw.h"
17#include "source_base/module_external/lapack_connector.h"
18#include "source_base/parallel_global.h"
19#include "source_io/module_parameter/parameter.h"
20#include "source_estate/elecstate.h"
21#include "source_basis/module_pw/pw_basis_k.h"
22#include "source_cell/module_symmetry/symmetry.h"
23#include "source_psi/psi_init.h"
24#include "source_pw/module_pwdft/structure_factor.h"
25#include "source_base/tool_title.h"
26#include "source_base/timer.h"
27
28#include <limits>
29
30template <typename T, typename Device>
31void Exx_Lip<T, Device>::cal_exx()
32{
33 ModuleBase::TITLE("Exx_Lip", "cal_exx");
34 ModuleBase::timer::tick("Exx_Lip", "cal_exx");
35
36 this->wf_wg_cal();
37 this->psi_cal();
38 for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik)
39 {
40 this->phi_cal(this->k_pack, ik);
41
42 this->judge_singularity(ik);
43 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) {
44 for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) {
45 this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] = T(0.0);
46 } }
47 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type)
48 {
49 this->sum2_factor = 0.0;
50 if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) {
51 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) {
52 for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) {
53 this->sum3[iw_l][iw_r] = T(0.0);
54 } } }
55 }
56
57 for (int iq_tmp = this->iq_vecik; iq_tmp < this->iq_vecik + this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin; ++iq_tmp) // !!! k_point
58 // parallel incompleted.need to loop iq in other pool
59 {
60 const int iq =
61 (ik < (this->k_pack->kv_ptr->get_nks() / PARAM.inp.nspin))
62 ? (iq_tmp % (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin))
63 : (iq_tmp % (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin) + (this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin));
64 this->qkg2_exp(ik, iq);
65 for (int ib = 0; ib < PARAM.inp.nbands; ++ib)
66 {
67 this->b_cal(ik, iq, ib);
68 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) {
69 if (iq == this->iq_vecik) {
70 this->sum3_cal(iq, ib);
71 } }
72 this->b_sum(iq, ib);
73 }
74 }
75 this->sum_all(ik);
76 }
77 this->exx_energy_cal();
78
79 ModuleBase::timer::tick("Exx_Lip", "cal_exx");
80}
81
82template <typename T, typename Device>
83Exx_Lip<T, Device>::Exx_Lip(const Exx_Info::Exx_Info_Lip& info_in,
84 const ModuleSymmetry::Symmetry& symm,
85 K_Vectors* kv_ptr_in,
86 psi::Psi<T, Device>* psi_local_in,
87 psi::Psi<T, Device>* kspw_psi_ptr_in,
88 const ModulePW::PW_Basis_K* wfc_basis_in,
89 const ModulePW::PW_Basis* rho_basis_in,
90 const Structure_Factor& sf,
91 const UnitCell* ucell_ptr_in,
92 const elecstate::ElecState* pelec_in)
93 : info(info_in)
94{
95 ModuleBase::TITLE("Exx_Lip", "init");
96 ModuleBase::timer::tick("Exx_Lip", "init");
97
98 this->k_pack = new k_package;
99 this->k_pack->kv_ptr = kv_ptr_in;
100 this->k_pack->psi_local = psi_local_in;
101 this->k_pack->pelec = pelec_in;
102 this->k_pack->kspw_psi_ptr = kspw_psi_ptr_in;
103 this->wfc_basis = wfc_basis_in;
104 this->rho_basis = rho_basis_in;
105 this->ucell_ptr = ucell_ptr_in;
106
107 int gzero_judge = -1;
108 if (this->rho_basis->gg_uniq[0] < 1e-8)
109 { gzero_judge = GlobalV::RANK_IN_POOL; }
110 #ifdef __MPI
111 MPI_Allreduce(&gzero_judge, &gzero_rank_in_pool, 1, MPI_INT, MPI_MAX, POOL_WORLD);
112 #endif
113 this->k_pack->wf_wg.create(this->k_pack->kv_ptr->get_nks(),PARAM.inp.nbands);
114
115 this->k_pack->hvec_array = new psi::Psi<T, Device>(this->k_pack->kv_ptr->get_nks(), PARAM.inp.nbands, PARAM.globalv.nlocal, kv_ptr_in->ngk, true);
116 // this->k_pack->hvec_array = new ModuleBase::ComplexMatrix[this->k_pack->kv_ptr->get_nks()];
117 // for( int ik=0; ik<this->k_pack->kv_ptr->get_nks(); ++ik)
118 // {
119 // this->k_pack->hvec_array[ik].create(PARAM.globalv.nlocal,PARAM.inp.nbands);
120 // }
121
122 // if (PARAM.inp.init_chg=="atomic")
123 {
124 this->q_pack = this->k_pack;
125 }
126 // else if(PARAM.inp.init_chg=="file")
127 // {
128 // read_q_pack(symm, this->wfc_basis, sf);
129 // }
130
131 this->phi.resize(PARAM.globalv.nlocal);
132 for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw)
133 { this->phi[iw].resize(this->rho_basis->nrxx); }
134
135 this->psi.resize(this->q_pack->kv_ptr->get_nks());
136 for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq)
137 {
138 this->psi[iq].resize(PARAM.inp.nbands);
139 for (int ib = 0; ib < PARAM.inp.nbands; ++ib)
140 { this->psi[iq][ib].resize(this->rho_basis->nrxx); }
141 }
142
143 this->recip_qkg2.resize(this->rho_basis->npw);
144
145 this->b.resize(PARAM.globalv.nlocal * this->rho_basis->npw);
146
147 this->sum1.resize(PARAM.globalv.nlocal * PARAM.globalv.nlocal);
148
149 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type)
150 {
151 if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)
152 {
153 this->b0.resize(PARAM.globalv.nlocal);
154 this->sum3.resize(PARAM.globalv.nlocal);
155 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l)
156 { this->sum3[iw_l].resize(PARAM.globalv.nlocal); }
157 }
158 }
159
160 this->exx_matrix.resize(this->k_pack->kv_ptr->get_nks());
161 for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik)
162 {
163 this->exx_matrix[ik].resize(PARAM.globalv.nlocal);
164 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l)
165 { this->exx_matrix[ik][iw_l].resize(PARAM.globalv.nlocal); }
166 }
167
168 ModuleBase::timer::tick("Exx_Lip", "init");
169}
170
171template <typename T, typename Device>
172Exx_Lip<T, Device>::~Exx_Lip()
173{
174 if (this->k_pack)
175 { delete this->k_pack->hvec_array; this->k_pack->hvec_array = nullptr; }
176 delete this->k_pack; this->k_pack = nullptr;
177
178 if (PARAM.inp.init_chg == "atomic")
179 {
180 this->q_pack = nullptr;
181 }
182 else if (PARAM.inp.init_chg == "file")
183 {
184 delete this->q_pack->kv_ptr; this->q_pack->kv_ptr = nullptr;
185 // delete this->q_pack->wf_ptr; this->q_pack->wf_ptr = nullptr;
186 // delete[] this->q_pack->hvec_array; this->q_pack->hvec_array=nullptr;
187 delete this->q_pack; this->q_pack = nullptr;
188 }
189}
190
191template <typename T, typename Device>
192void Exx_Lip<T, Device>::wf_wg_cal()
193{
194 ModuleBase::TITLE("Exx_Lip", "wf_wg_cal");
195 ModuleBase::timer::tick("Exx_Lip", "wf_wg_cal");
196 if (PARAM.inp.nspin == 1) {
197 for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) {
198 for (int ib = 0; ib < PARAM.inp.nbands; ++ib) {
199 this->k_pack->wf_wg(ik, ib) = this->k_pack->pelec->wg(ik, ib) / 2;
200 } } }
201 else if (PARAM.inp.nspin == 2) {
202 for (int ik = 0; ik < this->k_pack->kv_ptr->get_nks(); ++ik) {
203 for (int ib = 0; ib < PARAM.inp.nbands; ++ib) {
204 this->k_pack->wf_wg(ik, ib) = this->k_pack->pelec->wg(ik, ib);
205 } } }
206 ModuleBase::timer::tick("Exx_Lip", "wf_wg_cal");
207}
208
209template <typename T, typename Device>
210void Exx_Lip<T, Device>::phi_cal(k_package* kq_pack, const int ikq)
211{
212 ModuleBase::timer::tick("Exx_Lip", "phi_cal");
213 std::vector<T> porter (this->wfc_basis->nrxx);
214 for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw)
215 {
216 // this->wfc_basis->recip2real(&kq_pack->wf_ptr->wanf2[ikq](iw,0), porter.data(), ikq);
217 this->wfc_basis->recip2real(&(kq_pack->psi_local->operator()(ikq, iw, 0)), porter.data(), ikq);
218 int ir = 0;
219 for (int ix = 0; ix < this->rho_basis->nx; ++ix)
220 {
221 const Treal phase_x = kq_pack->kv_ptr->kvec_d[ikq].x * ix / this->rho_basis->nx;
222 for (int iy = 0; iy < this->rho_basis->ny; ++iy)
223 {
224 const Treal phase_xy = phase_x + kq_pack->kv_ptr->kvec_d[ikq].y * iy / this->rho_basis->ny;
225 for (int iz = this->rho_basis->startz_current; iz < this->rho_basis->startz_current + this->rho_basis->nplane; ++iz)
226 {
227 const Treal phase_xyz = phase_xy + kq_pack->kv_ptr->kvec_d[ikq].z * iz / this->rho_basis->nz;
228 const T exp_tmp = std::exp(phase_xyz * this->two_pi_i);
229 this->phi[iw][ir] = porter[ir] * exp_tmp;
230 ++ir;
231 }
232 }
233 }
234 }
235 ModuleBase::timer::tick("Exx_Lip", "phi_cal");
236}
237
238template <typename T, typename Device>
239void Exx_Lip<T, Device>::psi_cal()
240{
241 ModuleBase::TITLE("Exx_Lip", "psi_cal");
242 ModuleBase::timer::tick("Exx_Lip", "psi_cal");
243 if (PARAM.inp.init_chg == "atomic")
244 {
245 std::vector<T> porter (this->wfc_basis->nrxx);
246 for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq)
247 {
248 for (int ib = 0; ib < PARAM.inp.nbands; ++ib)
249 {
250 this->wfc_basis->recip2real(&(this->q_pack->kspw_psi_ptr->operator()(iq, ib, 0)), porter.data(), iq);
251
252 int ir = 0;
253 for (int ix = 0; ix < this->rho_basis->nx; ++ix)
254 {
255 const Treal phase_x = this->q_pack->kv_ptr->kvec_d[iq].x * ix / this->rho_basis->nx;
256 for (int iy = 0; iy < this->rho_basis->ny; ++iy)
257 {
258 const Treal phase_xy = phase_x + this->q_pack->kv_ptr->kvec_d[iq].y * iy / this->rho_basis->ny;
259 for (int iz = this->rho_basis->startz_current; iz < this->rho_basis->startz_current + this->rho_basis->nplane; ++iz)
260 {
261 const Treal phase_xyz = phase_xy + this->q_pack->kv_ptr->kvec_d[iq].z * iz / this->rho_basis->nz;
262 const T exp_tmp = std::exp(phase_xyz * this->two_pi_i);
263 this->psi[iq][ib][ir] = porter[ir] * exp_tmp;
264 ++ir;
265 }
266 }
267 }
268 }
269 }
270 }
271 else if (PARAM.inp.init_chg == "file")
272 {
273 for (int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq)
274 {
275 this->phi_cal(this->q_pack, iq);
276 for (int ib = 0; ib < PARAM.inp.nbands; ++ib)
277 {
278 ModuleBase::GlobalFunc::ZEROS(this->psi[iq][ib].data(), this->rho_basis->nrxx);
279 for (int iw = 0; iw < PARAM.globalv.nlocal; ++iw)
280 {
281 for (int ir = 0; ir < this->rho_basis->nrxx; ++ir)
282 {
283 this->psi[iq][ib][ir] += (*this->q_pack->hvec_array)(iq, ib, iw) * this->phi[iw][ir];
284 }
285 }
286 }
287 }
288 }
290 // !!! k_point parallel incompleted. need to loop iq in other pool
292 ModuleBase::timer::tick("Exx_Lip", "psi_cal");
293}
294
295template <typename T, typename Device>
296void Exx_Lip<T, Device>::judge_singularity(const int ik)
297{
298 ModuleBase::timer::tick("Exx_Lip", "judge_singularity");
299 if (PARAM.inp.init_chg=="atomic")
300 {
301 this->iq_vecik = ik;
302 }
303 else if(PARAM.inp.init_chg=="file")
304 {
305 Treal min_q_minus_k(std::numeric_limits<Treal>::max());
306 for( int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
307 {
308 const Treal q_minus_k((this->q_pack->kv_ptr->kvec_c[iq] - this->k_pack->kv_ptr->kvec_c[ik]).norm2());
309 if(q_minus_k < min_q_minus_k)
310 {
311 min_q_minus_k = q_minus_k;
312 this->iq_vecik = iq;
313 }
314 }
315 }
316 ModuleBase::timer::tick("Exx_Lip", "judge_singularity");
317}
318
319template <typename T, typename Device>
320void Exx_Lip<T, Device>::qkg2_exp(const int ik, const int iq)
321{
322 ModuleBase::timer::tick("Exx_Lip", "qkg2_exp");
323 for( int ig=0; ig<this->rho_basis->npw; ++ig)
324 {
325 const Treal qkg2 = ((this->q_pack->kv_ptr->kvec_c[iq] - this->k_pack->kv_ptr->kvec_c[ik] + this->rho_basis->gcar[ig]) * (ModuleBase::TWO_PI / this->ucell_ptr->lat0)).norm2();
326 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type)
327 {
328 if (std::abs(qkg2) < 1e-10)
329 { this->recip_qkg2[ig] = 0.0; } // 0 to ignore bb/qkg2 when qkg2==0
330 else
331 { this->recip_qkg2[ig] = 1.0 / qkg2; }
332 this->sum2_factor += this->recip_qkg2[ig] * std::exp(-info.lambda * qkg2);
333 this->recip_qkg2[ig] = sqrt(this->recip_qkg2[ig]);
334 }
335 else if (Conv_Coulomb_Pot_K::Ccp_Type::Erfc == info.ccp_type)
336 {
337 if (std::abs(qkg2) < 1e-10)
338 { this->recip_qkg2[ig] = 1.0 / (2 * info.hse_omega); }
339 else
340 { this->recip_qkg2[ig] = sqrt((1 - std::exp(-qkg2 / (4 * info.hse_omega * info.hse_omega))) / qkg2); }
341 }
342 else
343 {
344 throw( std::string(__FILE__) + " line " + std::to_string(__LINE__) );
345 }
346 }
347 ModuleBase::timer::tick("Exx_Lip", "qkg2_exp");
348}
349
350template <typename T, typename Device>
351void Exx_Lip<T, Device>::b_cal(const int ik, const int iq, const int ib)
352{
353 ModuleBase::timer::tick("Exx_Lip", "b_cal");
354 const ModuleBase::Vector3<double> q_minus_k = this->q_pack->kv_ptr->kvec_d[iq] - this->k_pack->kv_ptr->kvec_d[ik];
355 std::vector<T > mul_tmp(this->rho_basis->nrxx);
356 for( size_t ir=0,ix=0; ix<this->rho_basis->nx; ++ix)
357 {
358 const Treal phase_x = q_minus_k.x * ix / this->rho_basis->nx;
359 for( size_t iy=0; iy<this->rho_basis->ny; ++iy)
360 {
361 const Treal phase_xy = phase_x + q_minus_k.y * iy / this->rho_basis->ny;
362 for( size_t iz=this->rho_basis->startz_current; iz<this->rho_basis->startz_current+this->rho_basis->nplane; ++iz)
363 {
364 const Treal phase_xyz = phase_xy + q_minus_k.z * iz / this->rho_basis->nz;
365 mul_tmp[ir] = std::exp(-phase_xyz * this->two_pi_i);
366 mul_tmp[ir] *= this->psi[iq][ib][ir];
367 ++ir;
368 }
369 }
370 }
371
372 std::vector<T> porter (this->rho_basis->nrxx);
373 for(size_t iw=0; iw< PARAM.globalv.nlocal; ++iw)
374 {
375 auto& phi_w = this->phi[iw];
376 for( size_t ir=0; ir<this->rho_basis->nrxx; ++ir)
377 {
378 porter[ir] = conj(phi_w[ir]) * mul_tmp[ir] ;
379 // porter[ir] = phi_w[ir] * psi_q_b[ir] *exp_tmp[ir] ;
380 }
381 T* const b_w = &this->b[iw * this->rho_basis->npw];
382 this->rho_basis->real2recip( porter.data(), b_w);
383 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) {
384 if ((iq == this->iq_vecik) && (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)) {
385 this->b0[iw] = b_w[this->rho_basis->ig_gge0];
386 } }
387
388 for (size_t ig = 0; ig < this->rho_basis->npw; ++ig)
389 { b_w[ig] *= this->recip_qkg2[ig]; }
390 }
391 ModuleBase::timer::tick("Exx_Lip", "b_cal");
392}
393
394template <typename T, typename Device>
395void Exx_Lip<T, Device>::sum3_cal(const int iq, const int ib)
396{
397 ModuleBase::timer::tick("Exx_Lip", "sum3_cal");
398 if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) {
399 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l) {
400 for (int iw_r = 0; iw_r < PARAM.globalv.nlocal; ++iw_r) {
401 this->sum3[iw_l][iw_r] += this->b0[iw_l] * conj(this->b0[iw_r]) * (Treal)this->q_pack->wf_wg(iq, ib);
402 } } }
403 ModuleBase::timer::tick("Exx_Lip", "sum3_cal");
404}
405
406template <typename T, typename Device>
407void Exx_Lip<T, Device>::b_sum(const int iq, const int ib) // Peize Lin change 2019-04-14
408{
409 ModuleBase::timer::tick("Exx_Lip", "b_sum");
410 // this->sum1[iw_l,iw_r] += \sum_{ig} this->b[iw_l,ig] * conj(this->b[iw_r,ig]) * this->q_pack->wf_wg(iq,ib)
411 LapackConnector::herk(
412 'U','N',
413 PARAM.globalv.nlocal, this->rho_basis->npw,
414 (Treal)this->q_pack->wf_wg(iq, ib), this->b.data(), this->rho_basis->npw,
415 1.0, this->sum1.data(), PARAM.globalv.nlocal);
416 // cblas_zherk( CblasRowMajor, CblasUpper, CblasNoTrans,
417 // PARAM.globalv.nlocal, this->rho_basis->npw,
418 // this->q_pack->wf_wg(iq,ib), static_cast<void*>(this->b), this->rho_basis->npw,
419 // 1.0, static_cast<void*>(this->sum1), PARAM.globalv.nlocal);
420 ModuleBase::timer::tick("Exx_Lip", "b_sum");
421}
422
423template <typename T, typename Device>
424void Exx_Lip<T, Device>::sum_all(const int ik)
425{
426 ModuleBase::timer::tick("Exx_Lip", "sum_all");
427 Treal sum2_factor_g = 0.0;
428 const Treal fourpi_div_omega = 4 * (Treal)(ModuleBase::PI / this->ucell_ptr->omega);
429 const Treal spin_fac = 2.0;
430 #ifdef __MPI
431 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type)
432 { MPI_Reduce(&this->sum2_factor, &sum2_factor_g, 1, MPI_DOUBLE, MPI_SUM, gzero_rank_in_pool, POOL_WORLD); }
433 #endif
434 for (size_t iw_l = 1; iw_l < PARAM.globalv.nlocal; ++iw_l) {
435 for (size_t iw_r = 0; iw_r < iw_l; ++iw_r) {
436 this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] = conj(this->sum1[iw_r * PARAM.globalv.nlocal + iw_l]); // Peize Lin add conj 2019-04-14
437 } }
438
439 for (int iw_l = 0; iw_l < PARAM.globalv.nlocal; ++iw_l)
440 {
441 for( int iw_r=0; iw_r<PARAM.globalv.nlocal; ++iw_r)
442 {
443 this->exx_matrix[ik][iw_l][iw_r] = -fourpi_div_omega * this->sum1[iw_l * PARAM.globalv.nlocal + iw_r] * spin_fac;
444 if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type)
445 {
446 if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)
447 {
448 this->exx_matrix[ik][iw_l][iw_r] += spin_fac * (fourpi_div_omega * this->sum3[iw_l][iw_r] * sum2_factor_g);
449 this->exx_matrix[ik][iw_l][iw_r] += spin_fac * (-1 / (Treal)sqrt(info.lambda * ModuleBase::PI) * (Treal)(this->q_pack->kv_ptr->get_nks() / PARAM.inp.nspin) * this->sum3[iw_l][iw_r]);
450 }
451 }
452 }
453 }
454 ModuleBase::timer::tick("Exx_Lip", "sum_all");
455}
456
457template <typename T, typename Device>
458void Exx_Lip<T, Device>::exx_energy_cal()
459{
460 ModuleBase::TITLE("Exx_Lip","exx_energy_cal");
461 ModuleBase::timer::tick("Exx_Lip", "exx_energy_cal");
462
463 Treal exx_energy_tmp = 0.0;
464
465 for( int ik=0; ik<this->k_pack->kv_ptr->get_nks(); ++ik) {
466 for( int iw_l=0; iw_l<PARAM.globalv.nlocal; ++iw_l) {
467 for( int iw_r=0; iw_r<PARAM.globalv.nlocal; ++iw_r) {
468 for( int ib=0; ib<PARAM.inp.nbands; ++ib) {
469 exx_energy_tmp += (this->exx_matrix[ik][iw_l][iw_r] * conj((*this->k_pack->hvec_array)(ik, ib, iw_l)) * (*this->k_pack->hvec_array)(ik, ib, iw_r)).real() * this->k_pack->wf_wg(ik, ib);
470 } } } }
471 #ifdef __MPI
472 MPI_Allreduce( &exx_energy_tmp, &this->exx_energy, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); // !!! k_point parallel incompleted. different pools have different kv.set_nks(>) deadlock
473 #endif
474 this->exx_energy *= (PARAM.inp.nspin==1) ? 2 : 1;
475 this->exx_energy /= 2; // ETOT = E_band - 1/2 E_exx
476 ModuleBase::timer::tick("Exx_Lip", "exx_energy_cal");
477}
478
479template <typename T, typename Device>
480void Exx_Lip<T, Device>::write_q_pack() const
481{
482 ModuleBase::timer::tick("Exx_Lip", "write_q_pack");
483
484 if (PARAM.inp.out_chg[0] == 0)
485 { return; }
486
487 if (!GlobalV::RANK_IN_POOL)
488 {
489 const std::string exx_q_pack = "exx_q_pack/";
490
491 const std::string command_mkdir = "test -d " + PARAM.globalv.global_out_dir + exx_q_pack + " || mkdir " + PARAM.globalv.global_out_dir + exx_q_pack;
492 assert( system(command_mkdir.c_str()) == 0);
493
494 const std::string command_kpoint = "test -f " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file + " || cp " + PARAM.inp.kpoint_file + " " + PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file;
495 assert( system(command_kpoint.c_str()) == 0);
496
497 std::stringstream ss_wf_wg;
498 ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL;
499 std::ofstream ofs_wf_wg(ss_wf_wg.str().c_str());
500 for( int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq)
501 {
502 for( int ib=0; ib<PARAM.inp.nbands; ++ib)
503 {
504 ofs_wf_wg<<this->q_pack->wf_wg(iq,ib)<<"\t";
505 }
506 ofs_wf_wg<<std::endl;
507 }
508 ofs_wf_wg.close();
509
510 std::stringstream ss_hvec;
511 ss_hvec << PARAM.globalv.global_out_dir << exx_q_pack << "hvec_" << GlobalV::MY_POOL;
512 std::ofstream ofs_hvec(ss_hvec.str().c_str());
513 for( int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
514 {
515 for( int iw=0; iw<PARAM.globalv.nlocal; ++iw)
516 {
517 for( int ib=0; ib<PARAM.inp.nbands; ++ib)
518 {
519 ofs_hvec << (*this->q_pack->hvec_array)(iq, ib, iw).real() << " " << (*this->q_pack->hvec_array)(iq, ib, iw).imag() << " ";
520 }
521 ofs_hvec<<std::endl;
522 }
523 }
524 ofs_hvec.close();
525 }
526 ModuleBase::timer::tick("Exx_Lip", "write_q_pack");
527}
528
529/*
530void Exx_Lip::read_q_pack(const ModuleSymmetry::Symmetry& symm,
531 const ModulePW::PW_Basis_K* this->wfc_basis,
532 const Structure_Factor& sf)
533{
534 const std::string exx_q_pack = "exx_q_pack/";
535 this->q_pack = new k_package();
536 this->q_pack->kv_ptr = new K_Vectors();
537 const std::string exx_kpoint_card = PARAM.globalv.global_out_dir + exx_q_pack + PARAM.inp.kpoint_file;
538 this->q_pack->kv_ptr->set( symm, exx_kpoint_card, PARAM.inp.nspin, this->ucell_ptr->G, this->ucell_ptr->latvec, GlobalV::ofs_running );
539 this->q_pack->wf_ptr = new wavefunc();
540 this->q_pack->wf_ptr->allocate(this->q_pack->kv_ptr->get_nkstot(),
541 this->q_pack->kv_ptr->get_nks(),
542 this->q_pack->kv_ptr->ngk.data(),
543 this->wfc_basis->npwk_max); // mohan update 2021-02-25
544 // this->q_pack->wf_ptr->init(this->q_pack->kv_ptr->get_nks(),this->q_pack->kv_ptr,this->ucell_ptr,old_pwptr,&ppcell,&GlobalC::ORB,&hm,&Pkpoints);
545 this->q_pack->wf_ptr->table_local.create(ucell.ntype, ucell.nmax_total, PARAM.globalv.nqx);
546 // this->q_pack->wf_ptr->table_local.create(this->q_pack->wf_ptr->this->ucell_ptr->ntype, this->q_pack->wf_ptr->this->ucell_ptr->nmax_total, PARAM.globalv.nqx);
547 #ifdef __LCAO
548 Wavefunc_in_pw::make_table_q(GlobalC::ORB.orbital_file, this->q_pack->wf_ptr->table_local);
549 // Wavefunc_in_pw::make_table_q(this->q_pack->wf_ptr->ORB_ptr->orbital_file, this->q_pack->wf_ptr->table_local, this->q_pack->wf_ptr);
550 for(int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
551 {
552 Wavefunc_in_pw::produce_local_basis_in_pw(iq,
553 this->wfc_basis,
554 sf,
555 this->q_pack->wf_ptr->wanf2[iq],
556 this->q_pack->wf_ptr->table_local);
557 // Wavefunc_in_pw::produce_local_basis_in_pw(iq, this->q_pack->wf_ptr->wanf2[iq], this->q_pack->wf_ptr->table_local,
558 // this->q_pack->wf_ptr);
559 }
560 #endif
561 this->q_pack->wf_wg.create(this->q_pack->kv_ptr->get_nks(),PARAM.inp.nbands);
562 if(!GlobalV::RANK_IN_POOL)
563 {
564 std::stringstream ss_wf_wg;
565 ss_wf_wg << PARAM.globalv.global_out_dir << exx_q_pack << "wf_wg_" << GlobalV::MY_POOL;
566 std::ifstream ifs_wf_wg(ss_wf_wg.str().c_str());
567 for( int iq = 0; iq < this->q_pack->kv_ptr->get_nks(); ++iq)
568 {
569 for( int ib=0; ib<PARAM.inp.nbands; ++ib)
570 {
571 ifs_wf_wg>>this->q_pack->wf_wg(iq,ib);
572 }
573 }
574 ifs_wf_wg.close();
575 }
576 #ifdef __MPI
577 MPI_Bcast( this->q_pack->wf_wg.c, this->q_pack->kv_ptr->get_nks()*PARAM.inp.nbands, MPI_DOUBLE, 0, POOL_WORLD);
578 #endif
579 this->q_pack->hvec_array = new ModuleBase::ComplexMatrix [this->q_pack->kv_ptr->get_nks()];
580 for( int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
581 {
582 this->q_pack->hvec_array[iq].create(PARAM.globalv.nlocal,PARAM.inp.nbands);
583 }
584 if(!GlobalV::RANK_IN_POOL)
585 {
586 std::stringstream ss_hvec;
587 ss_hvec << PARAM.globalv.global_out_dir << exx_q_pack << "hvec_" << GlobalV::MY_POOL;
588 std::ifstream ifs_hvec(ss_hvec.str().c_str());
589 for( int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
590 {
591 for( int iw=0; iw<PARAM.globalv.nlocal; ++iw)
592 {
593 for( int ib=0; ib<PARAM.inp.nbands; ++ib)
594 {
595 double a,this->b;
596 ifs_hvec>>a>>this->b;
597 this->q_pack->hvec_array[iq](iw,ib) = {a,this->b};
598 }
599 }
600 }
601 ifs_hvec.close();
602 }
603 #ifdef __MPI
604 for( int iq=0; iq<this->q_pack->kv_ptr->get_nks(); ++iq)
605 {
606 MPI_Bcast( this->q_pack->hvec_array[iq].c, PARAM.globalv.nlocal*PARAM.inp.nbands, MPI_DOUBLE_COMPLEX, 0, POOL_WORLD);
607 }
608 #endif
609 return;
610}
611*/
612
613#endif
Exx_Lip::gzero_rank_in_pool
int gzero_rank_in_pool
Definition exx_lip.h:73
Exx_Lip::b
std::vector< T > b
Definition exx_lip.h:95
Exx_Lip::phi
std::vector< std::vector< T > > phi
Definition exx_lip.h:91
Exx_Lip::ucell_ptr
const UnitCell * ucell_ptr
Definition exx_lip.h:123
Exx_Lip::q_pack
struct Exx_Lip::k_package * q_pack
Exx_Lip::k_pack
struct Exx_Lip::k_package * k_pack
Exx_Lip::sum_all
void sum_all(const int ik)
Definition exx_lip.hpp:424
Exx_Lip::info
const Exx_Info::Exx_Info_Lip & info
Definition exx_lip.h:36
Exx_Lip::b_sum
void b_sum(const int iq, const int ib)
Definition exx_lip.hpp:407
Exx_Lip::cal_exx
void cal_exx()
Definition exx_lip.hpp:31
Exx_Lip::wfc_basis
const ModulePW::PW_Basis_K * wfc_basis
Definition exx_lip.h:121
Exx_Lip::write_q_pack
void write_q_pack() const
Definition exx_lip.hpp:480
Exx_Lip::exx_energy_cal
void exx_energy_cal()
Definition exx_lip.hpp:458
Exx_Lip::~Exx_Lip
~Exx_Lip()
Definition exx_lip.hpp:172
Exx_Lip::recip_qkg2
std::vector< Treal > recip_qkg2
Definition exx_lip.h:93
Exx_Lip::qkg2_exp
void qkg2_exp(const int ik, const int iq)
Definition exx_lip.hpp:320
Exx_Lip::Treal
typename GetTypeReal< T >::type Treal
Definition exx_lip.h:31
Exx_Lip::exx_matrix
std::vector< std::vector< std::vector< T > > > exx_matrix
Definition exx_lip.h:100
Exx_Lip::phi_cal
void phi_cal(k_package *kq_pack, const int ikq)
Definition exx_lip.hpp:210
Exx_Lip::sum1
std::vector< T > sum1
Definition exx_lip.h:97
Exx_Lip::wf_wg_cal
void wf_wg_cal()
Definition exx_lip.hpp:192
Exx_Lip::Exx_Lip
Exx_Lip(const Exx_Info::Exx_Info_Lip &info_in)
Exx_Lip::sum3
std::vector< std::vector< T > > sum3
Definition exx_lip.h:98
Exx_Lip::judge_singularity
void judge_singularity(const int ik)
Definition exx_lip.hpp:296
Exx_Lip::b0
std::vector< T > b0
Definition exx_lip.h:96
Exx_Lip::sum3_cal
void sum3_cal(const int iq, const int ib)
Definition exx_lip.hpp:395
Exx_Lip::b_cal
void b_cal(const int ik, int iq, const int ib)
Definition exx_lip.hpp:351
Exx_Lip::psi_cal
void psi_cal()
Definition exx_lip.hpp:239
Exx_Lip::rho_basis
const ModulePW::PW_Basis * rho_basis
Definition exx_lip.h:120
K_Vectors
Definition klist.h:13
K_Vectors::get_nks
int get_nks() const
Definition klist.h:68
K_Vectors::kvec_d
std::vector< ModuleBase::Vector3< double > > kvec_d
Cartesian coordinates of k points.
Definition klist.h:16
K_Vectors::ngk
std::vector< int > ngk
wk, weight of k points
Definition klist.h:20
LapackConnector::herk
static void herk(const char uplo, const char trans, const int n, const int k, const double alpha, const std::complex< double > *A, const int lda, const double beta, std::complex< double > *C, const int ldc)
Definition lapack_connector.h:453
ModuleBase::Vector3
3 elements vector
Definition vector3.h:22
ModuleBase::Vector3::x
T x
Definition vector3.h:24
ModuleBase::Vector3::y
T y
Definition vector3.h:25
ModuleBase::Vector3::z
T z
Definition vector3.h:26
ModuleBase::matrix::create
void create(const int nrow, const int ncol, const bool flag_zero=true)
Definition matrix.cpp:122
ModuleBase::timer::tick
static void tick(const std::string &class_name_in, const std::string &name_in)
Use twice at a time: the first time, set start_flag to false; the second time, calculate the time dur...
Definition timer.cpp:57
ModulePW::PW_Basis_K
Special pw_basis class. It includes different k-points.
Definition pw_basis_k.h:57
ModulePW::PW_Basis
A class which can convert a function of "r" to the corresponding linear superposition of plane waves ...
Definition pw_basis.h:56
ModulePW::PW_Basis::nrxx
int nrxx
Definition pw_basis.h:120
ModulePW::PW_Basis::npw
int npw
Definition pw_basis.h:115
ModulePW::PW_Basis::gg_uniq
double * gg_uniq
Definition pw_basis.h:160
ModuleSymmetry::Symmetry
Definition symmetry.h:16
Parameter::inp
const Input_para & inp
Definition parameter.h:26
Parameter::globalv
const System_para & globalv
Definition parameter.h:30
Structure_Factor
Definition structure_factor.h:11
UnitCell
Definition unitcell.h:16
elecstate::ElecState
Definition elecstate.h:15
psi::Psi
Definition psi.h:37
T
#define T
Definition exp.cpp:237
exx_lip.h
global_function.h
lapack_connector.h
Conv_Coulomb_Pot_K::Ccp_Type::Erfc
@ Erfc
Conv_Coulomb_Pot_K::Ccp_Type::Ccp
@ Ccp
Conv_Coulomb_Pot_K::Ccp_Type::Hf
@ Hf
GlobalV::MY_POOL
int MY_POOL
global index of pool (count in pool)
Definition global_variable.cpp:22
GlobalV::RANK_IN_POOL
int RANK_IN_POOL
global index of pool (count in process), my_rank in each pool
Definition global_variable.cpp:26
ModuleBase::GlobalFunc::ZEROS
void ZEROS(std::complex< T > *u, const TI n)
Definition global_function.h:109
ModuleBase::TWO_PI
const double TWO_PI
Definition constants.h:21
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
psi
Definition exx_lip.h:23
conj
double conj(double a)
Definition operator_lr_hxc.cpp:14
POOL_WORLD
MPI_Comm POOL_WORLD
Definition parallel_comm.cpp:6
parallel_global.h
PARAM
Parameter PARAM
Definition parameter.cpp:3
psi_init.h
pw_basis_k.h
klist.h
Exx_Info::Exx_Info_Lip
Definition exx_info.h:41
Exx_Info::Exx_Info_Lip::ccp_type
const Conv_Coulomb_Pot_K::Ccp_Type & ccp_type
Definition exx_info.h:42
Exx_Lip::k_package
Definition exx_lip.h:77
Exx_Lip::k_package::pelec
const elecstate::ElecState * pelec
Definition exx_lip.h:86
Exx_Lip::k_package::kspw_psi_ptr
psi::Psi< T, Device > * kspw_psi_ptr
PW wavefunction.
Definition exx_lip.h:80
Exx_Lip::k_package::kv_ptr
K_Vectors * kv_ptr
Definition exx_lip.h:78
Exx_Lip::k_package::psi_local
psi::Psi< T, Device > * psi_local
NAOs in PW.
Definition exx_lip.h:81
Exx_Lip::k_package::hvec_array
psi::Psi< T, Device > * hvec_array
LCAO wavefunction, the eigenvectors from lapack diagonalization.
Definition exx_lip.h:85
Exx_Lip::k_package::wf_wg
ModuleBase::matrix wf_wg
Definition exx_lip.h:82
Input_para::init_chg
std::string init_chg
"file","atomic"
Definition input_parameter.h:49
Input_para::out_chg
std::vector< int > out_chg
output charge density. 0: no; 1: yes
Definition input_parameter.h:366
Input_para::nspin
int nspin
LDA ; LSDA ; non-linear spin.
Definition input_parameter.h:84
Input_para::nbands
int nbands
number of bands
Definition input_parameter.h:75
Input_para::kpoint_file
std::string kpoint_file
file contains k-points – xiaohui modify 2015-02-01
Definition input_parameter.h:57
System_para::nlocal
int nlocal
total number of local basis.
Definition system_parameter.h:23
System_para::global_out_dir
std::string global_out_dir
global output directory
Definition system_parameter.h:42
structure_factor.h
tool_title.h
wavefunc_in_pw.h
Generated by doxygen 1.9.8