Skip to content

Commit 186c3a7

Browse files
mohanchenmohanchen
committed
merge newest version
2 parents 3c2f37e + aa52a32 commit 186c3a7

File tree

16 files changed

+343
-32
lines changed

16 files changed

+343
-32
lines changed

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
99
#include "module_hamilt_lcao/module_dftu/dftu.h"
1010
#include "module_io/berryphase.h"
11+
#include "module_io/cal_ldos.h"
1112
#include "module_io/cube_io.h"
1213
#include "module_io/dos_nao.h"
1314
#include "module_io/io_dmk.h"
@@ -419,6 +420,15 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
419420
this->p_hamilt);
420421
}
421422

423+
// out ldos
424+
if (PARAM.inp.out_ldos[0])
425+
{
426+
ModuleIO::Cal_ldos<TK>::cal_ldos_lcao(reinterpret_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec),
427+
this->psi[0],
428+
this->Pgrid,
429+
ucell);
430+
}
431+
422432
// 6) print out exchange-correlation potential
423433
if (PARAM.inp.out_mat_xc)
424434
{

‎source/module_esolver/esolver_ks_pw.cpp‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -841,13 +841,14 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
841841
GlobalV::ofs_running);
842842
}
843843

844-
// 2) out ldos
845-
if (PARAM.inp.out_ldos)
846-
{
847-
ModuleIO::cal_ldos(reinterpret_cast<elecstate::ElecStatePW<std::complex<double>>*>(this->pelec),
848-
this->psi[0],
849-
this->Pgrid,
850-
ucell);
844+
// out ldos
845+
if (PARAM.inp.out_ldos[0])
846+
{
847+
ModuleIO::Cal_ldos<std::complex<double>>::cal_ldos_pw(
848+
reinterpret_cast<elecstate::ElecStatePW<std::complex<double>>*>(this->pelec),
849+
this->psi[0],
850+
this->Pgrid,
851+
ucell);
851852
}
852853

853854
//! 3) Calculate the spillage value, used to generate numerical atomic orbitals

‎source/module_hamilt_lcao/module_deepks/deepks_vdpre.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -106,7 +106,7 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal,
106106
if (std::is_same<TK, std::complex<double>>::value)
107107
{
108108
const double arg
109-
= -(kvec_d[ik] * ModuleBase::Vector3<double>(dR1 - dR2)) * ModuleBase::TWO_PI;
109+
= (kvec_d[ik] * ModuleBase::Vector3<double>(dR1 - dR2)) * ModuleBase::TWO_PI;
110110
kphase = std::complex<double>(cos(arg), sin(arg));
111111
}
112112
TK_tensor* kpase_ptr = reinterpret_cast<TK_tensor*>(&kphase);
@@ -122,7 +122,7 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal,
122122
{
123123
TK_tensor tmp = overlap_1->get_value(iw1, ib + m1)
124124
* overlap_2->get_value(iw2, ib + m2) * *kpase_ptr;
125-
accessor[ik][iw1][iw2][inl][m1][m2] += tmp;
125+
accessor[ik][iw1_all][iw2_all][inl][m1][m2] += tmp;
126126
}
127127
}
128128
ib += nm;
@@ -282,7 +282,7 @@ void DeePKS_domain::prepare_phialpha(const int nlocal,
282282
for (int m1 = 0; m1 < nm; ++m1) // nm = 1 for s, 3 for p, 5 for d
283283
{
284284
TK_tensor tmp = overlap->get_value(iw1, ib + m1) * *kpase_ptr;
285-
accessor[iat][nl][ik][iw1_local][m1] += tmp;
285+
accessor[iat][nl][ik][iw1_all][m1] += tmp;
286286
}
287287
ib += nm;
288288
nl++;

‎source/module_io/cal_ldos.cpp‎

Lines changed: 99 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,17 @@
11
#include "cal_ldos.h"
22
#include "cube_io.h"
3+
#include "module_base/blas_connector.h"
4+
#include "module_base/scalapack_connector.h"
35

4-
void ModuleIO::cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
5-
const psi::Psi<std::complex<double>>& psi,
6-
const Parallel_Grid& pgrid,
7-
const UnitCell& ucell)
6+
#include <type_traits>
7+
8+
namespace ModuleIO
9+
{
10+
template <typename T>
11+
void Cal_ldos<T>::cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>* pelec,
12+
const psi::Psi<std::complex<double>>& psi,
13+
const Parallel_Grid& pgrid,
14+
const UnitCell& ucell)
815
{
916
// energy range for ldos (efermi as reference)
1017
const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0;
@@ -16,13 +23,13 @@ void ModuleIO::cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pele
1623
for (int ik = 0; ik < pelec->klist->get_nks(); ++ik)
1724
{
1825
psi.fix_k(ik);
19-
double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
26+
const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
2027
int nbands = psi.get_nbands();
2128

2229
for (int ib = 0; ib < nbands; ib++)
2330
{
2431
pelec->basis->recip2real(&psi(ib, 0), wfcr.data(), ik);
25-
double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
32+
const double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
2633
if (eigenval >= emin && eigenval <= emax)
2734
{
2835
for (int ir = 0; ir < pelec->basis->nrxx; ir++)
@@ -38,5 +45,90 @@ void ModuleIO::cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pele
3845
<< "LDOS_" << PARAM.inp.stm_bias << "eV"
3946
<< ".cube";
4047

41-
ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, 11, 0);
48+
const int precision = PARAM.inp.out_ldos[1];
49+
ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0);
50+
}
51+
52+
#ifdef __LCAO
53+
template <typename T>
54+
void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
55+
const psi::Psi<T>& psi,
56+
const Parallel_Grid& pgrid,
57+
const UnitCell& ucell)
58+
{
59+
// energy range for ldos (efermi as reference)
60+
const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0;
61+
const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0;
62+
63+
// calulate dm-like
64+
const int nbands_local = psi.get_nbands();
65+
const int nbasis_local = psi.get_nbasis();
66+
67+
// psi.T * wk * psi.conj()
68+
// result[ik](iw1,iw2) = \sum_{ib} psi[ik](ib,iw1).T * wk(k) * psi[ik](ib,iw2).conj()
69+
for (int ik = 0; ik < psi.get_nk(); ++ik)
70+
{
71+
psi.fix_k(ik);
72+
const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
73+
74+
// T* dmk_pointer = DM.get_DMK_pointer(ik);
75+
76+
psi::Psi<T> wk_psi(1, psi.get_nbands(), psi.get_nbasis(), psi.get_nbasis(), true);
77+
const T* ppsi = psi.get_pointer();
78+
T* pwk_psi = wk_psi.get_pointer();
79+
80+
// #ifdef _OPENMP
81+
// #pragma omp parallel for schedule(static, 1024)
82+
// #endif
83+
// for (int i = 0; i < wk_psi.size(); ++i)
84+
// {
85+
// pwk_psi[i] = my_conj(ppsi[i]);
86+
// }
87+
88+
// int ib_global = 0;
89+
// for (int ib_local = 0; ib_local < nbands_local; ++ib_local)
90+
// {
91+
// while (ib_local != ParaV->global2local_col(ib_global))
92+
// {
93+
// ++ib_global;
94+
// if (ib_global >= wg.nc)
95+
// {
96+
// ModuleBase::WARNING_QUIT("cal_ldos", "please check global2local_col!");
97+
// }
98+
// }
99+
100+
// const double eigenval = (pelec->ekb(ik, ib_global) - efermi) * ModuleBase::Ry_to_eV;
101+
// if (eigenval >= emin && eigenval <= emax)
102+
// {
103+
// for (int ir = 0; ir < pelec->basis->nrxx; ir++)
104+
// ldos[ir] += pelec->klist->wk[ik] * norm(wfcr[ir]);
105+
// }
106+
107+
// double* wg_wfc_pointer = &(wk_psi(0, ib_local, 0));
108+
// BlasConnector::scal(nbasis_local, pelec->klist->wk[ik], wg_wfc_pointer, 1);
109+
// }
110+
111+
// // C++: dm(iw1,iw2) = psi(ib,iw1).T * wk_psi(ib,iw2)
112+
// #ifdef __MPI
113+
// psiMulPsiMpi(wk_psi, psi, dmk_pointer, ParaV->desc_wfc, ParaV->desc);
114+
// #else
115+
// psiMulPsi(wk_psi, psi, dmk_pointer);
116+
// #endif
117+
}
118+
}
119+
120+
double my_conj(double x)
121+
{
122+
return x;
42123
}
124+
125+
std::complex<double> my_conj(const std::complex<double>& z)
126+
{
127+
return {z.real(), -z.imag()};
128+
}
129+
130+
#endif
131+
132+
template class Cal_ldos<double>; // Gamma_only case
133+
template class Cal_ldos<std::complex<double>>; // multi-k case
134+
} // namespace ModuleIO

‎source/module_io/cal_ldos.h‎

Lines changed: 17 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,29 @@
11
#ifndef CAL_LDOS_H
22
#define CAL_LDOS_H
33

4+
#include "module_elecstate/elecstate_lcao.h"
45
#include "module_elecstate/elecstate_pw.h"
56

67
namespace ModuleIO
78
{
9+
template <typename T>
10+
class Cal_ldos
11+
{
12+
public:
13+
Cal_ldos(){};
14+
~Cal_ldos(){};
15+
16+
static void cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>* pelec,
17+
const psi::Psi<std::complex<double>>& psi,
18+
const Parallel_Grid& pgrid,
19+
const UnitCell& ucell);
820

9-
void cal_ldos(const elecstate::ElecStatePW<std::complex<double>>* pelec,
10-
const psi::Psi<std::complex<double>>& psi,
11-
const Parallel_Grid& pgrid,
12-
const UnitCell& ucell);
21+
static void cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
22+
const psi::Psi<T>& psi,
23+
const Parallel_Grid& pgrid,
24+
const UnitCell& ucell);
1325

26+
}; // namespace Cal_ldos
1427
} // namespace ModuleIO
1528

1629
#endif // CAL_LDOS_H

‎source/module_io/read_input_item_output.cpp‎

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -145,8 +145,17 @@ void ReadInput::item_output()
145145
}
146146
{
147147
Input_Item item("out_ldos");
148-
item.annotation = "output local density of states";
149-
read_sync_bool(input.out_ldos);
148+
item.annotation = "output local density of states, second parameter controls the precision";
149+
item.read_value = [](const Input_Item& item, Parameter& para) {
150+
const size_t count = item.get_size();
151+
if (count != 1 && count != 2)
152+
{
153+
ModuleBase::WARNING_QUIT("ReadInput", "out_ldos should have 1 or 2 values");
154+
}
155+
para.input.out_ldos[0] = assume_as_boolean(item.str_values[0]);
156+
para.input.out_ldos[1] = (count == 2) ? std::stoi(item.str_values[1]) : 3;
157+
};
158+
sync_intvec(input.out_ldos, 2, 0);
150159
this->add_item(item);
151160
}
152161
{

‎source/module_io/read_input_item_postprocess.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ void ReadInput::item_postprocess()
5858
item.annotation = "bias voltage used to calculate ldos";
5959
read_sync_double(input.stm_bias);
6060
item.check_value = [](const Input_Item& item, const Parameter& para) {
61-
if (para.input.out_ldos && para.input.stm_bias == 0.0)
61+
if (para.input.out_ldos[0] && para.input.stm_bias == 0.0)
6262
{
6363
ModuleBase::WARNING_QUIT("ReadInput", "a nonzero stm_bias is required for ldos calculation");
6464
}

‎source/module_io/test/read_input_ptest.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -201,7 +201,8 @@ TEST_F(InputParaTest, ParaRead)
201201
EXPECT_EQ(param.inp.out_wfc_pw, 0);
202202
EXPECT_EQ(param.inp.out_wfc_r, 0);
203203
EXPECT_EQ(param.inp.out_dos, 0);
204-
EXPECT_EQ(param.inp.out_ldos, true);
204+
EXPECT_EQ(param.inp.out_ldos[0], 1);
205+
EXPECT_EQ(param.inp.out_ldos[1], 3);
205206
EXPECT_EQ(param.inp.out_band[0], 0);
206207
EXPECT_EQ(param.inp.out_band[1], 8);
207208
EXPECT_EQ(param.inp.out_proj_band, 0);

‎source/module_io/test/support/INPUT‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,7 @@ out_pot 2 #output realspace potential
6565
out_wfc_pw 0 #output wave functions
6666
out_wfc_r 0 #output wave functions in realspace
6767
out_dos 0 #output energy and dos
68-
out_ldos True #output local density of states
68+
out_ldos 1 #output local density of states, second parameter controls the precision
6969
out_band 0 #output energy and band structure
7070
out_proj_band FaLse #output projected band structure
7171
restart_save f #print to disk every step for restart

‎source/module_io/test_serial/read_input_item_test.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -393,7 +393,7 @@ TEST_F(InputTest, Item_test)
393393
}
394394
{ // stm_bias
395395
auto it = find_label("stm_bias", readinput.input_lists);
396-
param.input.out_ldos = true;
396+
param.input.out_ldos[0] = 1;
397397
param.input.stm_bias = 0.0;
398398

399399
testing::internal::CaptureStdout();

0 commit comments

Comments
 (0)