Skip to content

Commit 48434f0

Browse files
YuLiu98YuLiu98
authored
Feature: ldos in lcao basis sets (#6123)
* Feature: ldos in lcao basis sets * Feature: enable multiple bias value in a scf * update unittests
1 parent ef95b28 commit 48434f0

File tree

7 files changed

+136
-103
lines changed

7 files changed

+136
-103
lines changed

‎docs/advanced/input_files/input-main.md‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1976,8 +1976,12 @@ These variables are used to control the calculation of DOS. [Detailed introducti
19761976

19771977
### stm_bias
19781978

1979-
- **Type**: Real
1980-
- **Description**: The bias voltage used to calculate local density of states to simulate scanning tunneling microscope, see details in [out_ldos](#out_ldos).
1979+
- **Type**: Real Real(optional) Integer(optional)
1980+
- **Description**: The bias voltage used to calculate local density of states to simulate scanning tunneling microscope, see details in [out_ldos](#out_ldos). When using three parameters:
1981+
1982+
- The first parameter specifies the initial bias voltage value.
1983+
- The second parameter defines the voltage increment (step size between consecutive bias values).
1984+
- The third parameter determines the total number of voltage points
19811985
- **Default**: 1.0
19821986
- **Unit**: V
19831987

‎source/module_io/cal_ldos.cpp‎

Lines changed: 105 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,8 @@
11
#include "cal_ldos.h"
2+
23
#include "cube_io.h"
3-
#include "module_base/blas_connector.h"
4-
#include "module_base/scalapack_connector.h"
4+
#include "module_elecstate/module_dm/cal_dm_psi.h"
5+
#include "module_hamilt_lcao/module_gint/temp_gint/gint_interface.h"
56

67
#include <type_traits>
78

@@ -13,40 +14,43 @@ void Cal_ldos<T>::cal_ldos_pw(const elecstate::ElecStatePW<std::complex<double>>
1314
const Parallel_Grid& pgrid,
1415
const UnitCell& ucell)
1516
{
16-
// energy range for ldos (efermi as reference)
17-
const double emin = PARAM.inp.stm_bias < 0 ? PARAM.inp.stm_bias : 0;
18-
const double emax = PARAM.inp.stm_bias > 0 ? PARAM.inp.stm_bias : 0;
19-
20-
std::vector<double> ldos(pelec->charge->nrxx);
21-
std::vector<std::complex<double>> wfcr(pelec->basis->nrxx);
22-
23-
for (int ik = 0; ik < pelec->klist->get_nks(); ++ik)
17+
for (int ie = 0; ie < PARAM.inp.stm_bias[2]; ie++)
2418
{
25-
psi.fix_k(ik);
26-
const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
27-
int nbands = psi.get_nbands();
19+
// energy range for ldos (efermi as reference)
20+
const double en = PARAM.inp.stm_bias[0] + ie * PARAM.inp.stm_bias[1];
21+
const double emin = en < 0 ? en : 0;
22+
const double emax = en > 0 ? en : 0;
23+
24+
std::vector<double> ldos(pelec->charge->nrxx);
25+
std::vector<std::complex<double>> wfcr(pelec->basis->nrxx);
2826

29-
for (int ib = 0; ib < nbands; ib++)
27+
for (int ik = 0; ik < pelec->klist->get_nks(); ++ik)
3028
{
31-
pelec->basis->recip2real(&psi(ib, 0), wfcr.data(), ik);
32-
const double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
33-
if (eigenval >= emin && eigenval <= emax)
29+
psi.fix_k(ik);
30+
const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
31+
int nbands = psi.get_nbands();
32+
33+
for (int ib = 0; ib < nbands; ib++)
3434
{
35-
for (int ir = 0; ir < pelec->basis->nrxx; ir++)
36-
{
37-
ldos[ir] += pelec->klist->wk[ik] * norm(wfcr[ir]);
38-
}
35+
pelec->basis->recip2real(&psi(ib, 0), wfcr.data(), ik);
36+
const double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
37+
if (eigenval >= emin && eigenval <= emax)
38+
{
39+
for (int ir = 0; ir < pelec->basis->nrxx; ir++)
40+
{
41+
ldos[ir] += pelec->klist->wk[ik] * norm(wfcr[ir]);
42+
}
43+
}
3944
}
4045
}
41-
}
4246

43-
std::stringstream fn;
44-
fn << PARAM.globalv.global_out_dir
45-
<< "LDOS_" << PARAM.inp.stm_bias << "eV"
46-
<< ".cube";
47+
std::stringstream fn;
48+
fn << PARAM.globalv.global_out_dir << "LDOS_" << en << "eV"
49+
<< ".cube";
4750

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);
51+
const int precision = PARAM.inp.out_ldos[1];
52+
ModuleIO::write_vdata_palgrid(pgrid, ldos.data(), 0, PARAM.inp.nspin, 0, fn.str(), 0, &ucell, precision, 0);
53+
}
5054
}
5155

5256
#ifdef __LCAO
@@ -56,75 +60,83 @@ void Cal_ldos<T>::cal_ldos_lcao(const elecstate::ElecStateLCAO<T>* pelec,
5660
const Parallel_Grid& pgrid,
5761
const UnitCell& ucell)
5862
{
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;
63+
for (int ie = 0; ie < PARAM.inp.stm_bias[2]; ie++)
64+
{
65+
// energy range for ldos (efermi as reference)
66+
const double en = PARAM.inp.stm_bias[0] + ie * PARAM.inp.stm_bias[1];
67+
const double emin = en < 0 ? en : 0;
68+
const double emax = en > 0 ? en : 0;
69+
70+
// calculate weight (for bands not in the range, weight is zero)
71+
ModuleBase::matrix weight(pelec->ekb.nr, pelec->ekb.nc);
72+
for (int ik = 0; ik < pelec->ekb.nr; ++ik)
73+
{
74+
const double efermi = pelec->eferm.get_efval(pelec->klist->isk[ik]);
6275

63-
// calulate dm-like
64-
const int nbands_local = psi.get_nbands();
65-
const int nbasis_local = psi.get_nbasis();
76+
for (int ib = 0; ib < pelec->ekb.nc; ib++)
77+
{
78+
const double eigenval = (pelec->ekb(ik, ib) - efermi) * ModuleBase::Ry_to_eV;
79+
if (eigenval >= emin && eigenval <= emax)
80+
{
81+
weight(ik, ib) = pelec->klist->wk[ik];
82+
}
83+
}
84+
}
6685

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-
}
86+
// calculate dm-like for ldos
87+
const int nspin_dm = PARAM.inp.nspin == 2 ? 2 : 1;
88+
elecstate::DensityMatrix<T, double> dm_ldos(pelec->DM->get_paraV_pointer(),
89+
nspin_dm,
90+
pelec->klist->kvec_d,
91+
pelec->klist->get_nks() / nspin_dm);
92+
93+
elecstate::cal_dm_psi(pelec->DM->get_paraV_pointer(), weight, psi, dm_ldos);
94+
dm_ldos.init_DMR(*(pelec->DM->get_DMR_pointer(1)));
95+
dm_ldos.cal_DMR();
96+
97+
// allocate ldos space
98+
std::vector<double> ldos_space(PARAM.inp.nspin * pelec->charge->nrxx);
99+
double** ldos = new double*[PARAM.inp.nspin];
100+
for (int is = 0; is < PARAM.inp.nspin; ++is)
101+
{
102+
ldos[is] = &ldos_space[is * pelec->charge->nrxx];
103+
}
119104

120-
double my_conj(double x)
121-
{
122-
return x;
123-
}
105+
// calculate ldos
106+
#ifndef __NEW_GINT
107+
ModuleBase::WARNING_QUIT("Cal_ldos::dm2ldos",
108+
"do not support old grid integral, please recompile with __NEW_GINT");
109+
#else
110+
ModuleGint::cal_gint_rho(dm_ldos.get_DMR_vector(), PARAM.inp.nspin, ldos);
111+
#endif
124112

125-
std::complex<double> my_conj(const std::complex<double>& z)
126-
{
127-
return {z.real(), -z.imag()};
113+
// I'm not sure whether ldos should be output for each spin or not
114+
// ldos[0] += ldos[1] for nspin_dm == 2
115+
if (nspin_dm == 2)
116+
{
117+
BlasConnector::axpy(pelec->charge->nrxx, 1.0, ldos[1], 1, ldos[0], 1);
118+
}
119+
120+
// write ldos to cube file
121+
std::stringstream fn;
122+
fn << PARAM.globalv.global_out_dir << "LDOS_" << en << "eV"
123+
<< ".cube";
124+
125+
const int precision = PARAM.inp.out_ldos[1];
126+
ModuleIO::write_vdata_palgrid(pgrid,
127+
ldos_space.data(),
128+
0,
129+
PARAM.inp.nspin,
130+
0,
131+
fn.str(),
132+
0,
133+
&ucell,
134+
precision,
135+
0);
136+
137+
// free memory
138+
delete[] ldos;
139+
}
128140
}
129141

130142
#endif

‎source/module_io/cal_ldos.h‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,6 @@ class Cal_ldos
2222
const psi::Psi<T>& psi,
2323
const Parallel_Grid& pgrid,
2424
const UnitCell& ucell);
25-
2625
}; // namespace Cal_ldos
2726
} // namespace ModuleIO
2827

‎source/module_io/read_input_item_postprocess.cpp‎

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -56,13 +56,27 @@ void ReadInput::item_postprocess()
5656
{
5757
Input_Item item("stm_bias");
5858
item.annotation = "bias voltage used to calculate ldos";
59-
read_sync_double(input.stm_bias);
59+
item.read_value = [](const Input_Item& item, Parameter& para) {
60+
const size_t count = item.get_size();
61+
if (count != 1 && count != 3)
62+
{
63+
ModuleBase::WARNING_QUIT("ReadInput", "stm_bias should have 1 or 3 values");
64+
}
65+
para.input.stm_bias[0] = std::stod(item.str_values[0]);
66+
para.input.stm_bias[1] = (count == 3) ? std::stod(item.str_values[1]) : 0.1;
67+
para.input.stm_bias[2] = (count == 3) ? std::stod(item.str_values[2]) : 1;
68+
};
6069
item.check_value = [](const Input_Item& item, const Parameter& para) {
61-
if (para.input.out_ldos[0] && para.input.stm_bias == 0.0)
70+
if (para.input.stm_bias[2] <= 0)
71+
{
72+
ModuleBase::WARNING_QUIT("ReadInput", "stm_bias[2] should be greater than 0");
73+
}
74+
if (para.input.stm_bias[1] == 0)
6275
{
63-
ModuleBase::WARNING_QUIT("ReadInput", "a nonzero stm_bias is required for ldos calculation");
76+
ModuleBase::WARNING_QUIT("ReadInput", "stm_bias[1] should be nonzero");
6477
}
6578
};
79+
sync_doublevec(input.stm_bias, 3, 0);
6680
this->add_item(item);
6781
}
6882

‎source/module_io/test/read_input_ptest.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -222,7 +222,7 @@ TEST_F(InputParaTest, ParaRead)
222222
EXPECT_DOUBLE_EQ(param.inp.dos_edelta_ev, 0.01);
223223
EXPECT_DOUBLE_EQ(param.inp.dos_scale, 0.01);
224224
EXPECT_DOUBLE_EQ(param.inp.dos_sigma, 0.07);
225-
EXPECT_DOUBLE_EQ(param.inp.stm_bias, 2.0);
225+
EXPECT_DOUBLE_EQ(param.inp.stm_bias[0], 2.0);
226226
EXPECT_FALSE(param.inp.out_element_info);
227227
EXPECT_DOUBLE_EQ(param.inp.lcao_ecut, 20);
228228
EXPECT_DOUBLE_EQ(param.inp.lcao_dk, 0.01);

‎source/module_io/test_serial/read_input_item_test.cpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -393,9 +393,13 @@ TEST_F(InputTest, Item_test)
393393
}
394394
{ // stm_bias
395395
auto it = find_label("stm_bias", readinput.input_lists);
396-
param.input.out_ldos[0] = 1;
397-
param.input.stm_bias = 0.0;
396+
param.input.stm_bias = {1.0, 0.1, 0.0};
397+
testing::internal::CaptureStdout();
398+
EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
399+
output = testing::internal::GetCapturedStdout();
400+
EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
398401

402+
param.input.stm_bias = {1.0, 0.0, 2.0};
399403
testing::internal::CaptureStdout();
400404
EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
401405
output = testing::internal::GetCapturedStdout();

‎source/module_parameter/input_parameter.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -418,7 +418,7 @@ struct Input_para
418418
double dos_scale = 0.01;
419419
double dos_sigma = 0.07; ///< pengfei 2014-10-13
420420
int dos_nche = 100; ///< orders of Chebyshev expansions for dos
421-
double stm_bias = 1.0; ///< bias voltage for STM
421+
std::vector<double> stm_bias = {1.0, 0.1, 1}; ///< bias voltage for STM (start value, step, number)
422422

423423
bool cal_cond = false; ///< calculate electronic conductivities
424424
double cond_che_thr = 1e-8; ///< control the error of Chebyshev expansions

0 commit comments

Comments
 (0)