Skip to content

Commit edb2d1c

Browse files
ESROAMERmohanchen
ESROAMER
and
mohanchen
authored
Fix HR output in tddft (#4150)
* outHR * Delete source/LM * Delete abacus_cmake.sh * Update H_TDDFT_pw.h * Update td_ekinetic_lcao.cpp * update HR output * Delete abacus_cmake.sh * Delete source/LM * Update H_TDDFT_pw.h * Update LCAO_matrix.cpp * Update Makefile.Objects * Update esolver_ks_lcao_tddft.cpp * update * Delete abacus_cmake.sh * fix out_mat_R * Update for_testing_input_conv.h * remove GlobalC and GlobalV int td operators --------- Co-authored-by: Mohan Chen <[email protected]>
1 parent 840f0d5 commit edb2d1c

File tree

16 files changed

+255
-34
lines changed

16 files changed

+255
-34
lines changed

‎source/Makefile.Objects

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -498,6 +498,7 @@ OBJS_LCAO=DM_gamma.o\
498498
middle_hamilt.o\
499499
norm_psi.o\
500500
propagator.o\
501+
td_velocity.o\
501502
upsi.o\
502503
FORCE_STRESS.o\
503504
FORCE_gamma.o\

‎source/module_elecstate/potentials/H_TDDFT_pw.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -133,4 +133,4 @@ class H_TDDFT_pw : public PotBase
133133

134134
} // namespace elecstate
135135

136-
#endif
136+
#endif

‎source/module_esolver/esolver_ks_lcao_tddft.cpp

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include "module_elecstate/module_charge/symmetry_rho.h"
1616
#include "module_elecstate/occupy.h"
1717
#include "module_hamilt_lcao/module_tddft/evolve_elec.h"
18+
#include "module_hamilt_lcao/module_tddft/td_velocity.h"
1819
#include "module_hamilt_pw/hamilt_pwdft/global.h"
1920
#include "module_io/print_info.h"
2021

‎source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp

Lines changed: 47 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include "module_hamilt_lcao/hamilt_lcaodft/center2_orb-orb11.h"
66
#include "module_elecstate/potentials/H_TDDFT_pw.h"
77
#include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h"
8+
#include "module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h"
89

910
#include "module_cell/module_neighbor/sltk_grid_driver.h"
1011
#include "module_base/libm/libm.h"
@@ -29,6 +30,10 @@ TDEkinetic<OperatorLCAO<TK, TR>>::TDEkinetic(LCAO_Matrix* LM_in,
2930
this->init_td();
3031
// initialize HR to get adjs info.
3132
this->initialize_HR(Grid,this->LM->ParaV);
33+
if(TD_Velocity::out_mat_R == true)
34+
{
35+
out_mat_R = true;
36+
}
3237
}
3338
template <typename TK, typename TR>
3439
TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
@@ -37,6 +42,7 @@ TDEkinetic<OperatorLCAO<TK, TR>>::~TDEkinetic()
3742
{
3843
delete this->hR_tmp;
3944
}
45+
TD_Velocity::td_vel_op = nullptr;
4046
}
4147
//term A^2*S
4248
template <typename TK, typename TR>
@@ -185,22 +191,24 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::cal_HR_IJR(const int& iat1,
185191
template <typename TK, typename TR>
186192
void TDEkinetic<OperatorLCAO<TK, TR>>::init_td(void)
187193
{
194+
TD_Velocity::td_vel_op = &td_velocity;
188195
//calculate At in cartesian coorinates.
189-
double l_norm[3]={GlobalC::ucell.a1.norm() ,GlobalC::ucell.a2.norm() ,GlobalC::ucell.a3.norm()};
196+
double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
190197
double (&A)[3] = elecstate::H_TDDFT_pw::At;
191-
cart_At = GlobalC::ucell.a1*A[0]/l_norm[0] + GlobalC::ucell.a2*A[1]/l_norm[1] + GlobalC::ucell.a3*A[2]/l_norm[2];
198+
cart_At = this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2];
192199
std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
193-
200+
194201
//init MOT,MGT
202+
const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance();
195203
this->MOT.allocate(
196-
GlobalC::ORB.get_ntype(), // number of atom types
197-
GlobalC::ORB.get_lmax(), // max L used to calculate overlap
198-
static_cast<int>(GlobalC::ORB.get_kmesh()) | 1, // kpoints, for integration in k space
199-
GlobalC::ORB.get_Rmax(), // max value of radial table
200-
GlobalC::ORB.get_dR(), // delta R, for making radial table
201-
GlobalC::ORB.get_dk()); // Peize Lin change 2017-04-16
204+
orb.get_ntype(), // number of atom types
205+
orb.get_lmax(), // max L used to calculate overlap
206+
static_cast<int>(orb.get_kmesh()) | 1, // kpoints, for integration in k space
207+
orb.get_Rmax(), // max value of radial table
208+
orb.get_dR(), // delta R, for making radial table
209+
orb.get_dk()); // Peize Lin change 2017-04-16
202210
int Lmax_used, Lmax;
203-
this->MOT.init_Table_Spherical_Bessel (2, 1, Lmax_used, Lmax, 1, GlobalC::ORB, GlobalC::ucell.infoNL.Beta);
211+
this->MOT.init_Table_Spherical_Bessel (2, 1, Lmax_used, Lmax, 1, orb, this->ucell->infoNL.Beta);
204212

205213
//=========================================
206214
// (2) init Ylm Coef
@@ -214,17 +222,17 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::init_td(void)
214222
this->MGT.init_Gaunt( Lmax );
215223

216224
//init_radial table
217-
for( size_t TA=0; TA!=GlobalC::ORB.get_ntype(); ++TA )
218-
for( size_t TB=0; TB!=GlobalC::ORB.get_ntype(); ++TB )
219-
for( int LA=0; LA<=GlobalC::ORB.Phi[TA].getLmax(); ++LA )
220-
for( size_t NA=0; NA!=GlobalC::ORB.Phi[TA].getNchi(LA); ++NA )
221-
for( int LB=0; LB<=GlobalC::ORB.Phi[TB].getLmax(); ++LB )
222-
for( size_t NB=0; NB!=GlobalC::ORB.Phi[TB].getNchi(LB); ++NB )
225+
for( size_t TA=0; TA!=orb.get_ntype(); ++TA )
226+
for( size_t TB=0; TB!=orb.get_ntype(); ++TB )
227+
for( int LA=0; LA<=orb.Phi[TA].getLmax(); ++LA )
228+
for( size_t NA=0; NA!=orb.Phi[TA].getNchi(LA); ++NA )
229+
for( int LB=0; LB<=orb.Phi[TB].getLmax(); ++LB )
230+
for( size_t NB=0; NB!=orb.Phi[TB].getNchi(LB); ++NB )
223231
center2_orb11_s[TA][TB][LA][NA][LB].insert(
224232
std::make_pair(NB,
225233
Center2_Orb::Orb11(
226-
GlobalC::ORB.Phi[TA].PhiLN(LA,NA),
227-
GlobalC::ORB.Phi[TB].PhiLN(LB,NB),
234+
orb.Phi[TA].PhiLN(LA,NA),
235+
orb.Phi[TB].PhiLN(LB,NB),
228236
this->MOT, this->MGT)));
229237
for( auto &coA : center2_orb11_s )
230238
for( auto &coB : coA.second )
@@ -361,22 +369,39 @@ void TDEkinetic<OperatorLCAO<TK, TR>>::contributeHk(int ik)
361369
template<>
362370
void TDEkinetic<OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik)
363371
{
364-
if (GlobalV::ESOLVER_TYPE != "tddft" || elecstate::H_TDDFT_pw::stype != 1)
372+
if (TD_Velocity::tddft_velocity == false)
365373
{
366374
return;
367375
}
368376
else{
369377
ModuleBase::TITLE("TDEkinetic", "contributeHk");
370378
ModuleBase::timer::tick("TDEkinetic", "contributeHk");
379+
const Parallel_Orbitals* paraV = this->hR_tmp->get_atom_pair(0).get_paraV();
380+
//save HR data for output
381+
int spin_tot = paraV->nspin;
382+
if(spin_tot==4);
383+
else if(!output_hR_done && out_mat_R)
384+
{
385+
for(int spin_now = 0;spin_now < spin_tot;spin_now++)
386+
{
387+
sparse_format::cal_HContainer_cd(
388+
*(paraV),
389+
spin_now,
390+
1e-10,
391+
*hR_tmp,
392+
td_velocity.HR_sparse_td_vel[spin_now]);
393+
}
394+
output_hR_done = true;
395+
}
371396
//folding inside HR to HK
372397
if(ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER())
373398
{
374-
const int nrow = this->LM->ParaV->get_row_size();
399+
const int nrow = paraV->get_row_size();
375400
hamilt::folding_HR(*this->hR_tmp, this->hK->data(), this->kvec_d[ik], nrow, 1);
376401
}
377402
else
378403
{
379-
const int ncol = this->LM->ParaV->get_col_size();
404+
const int ncol = paraV->get_col_size();
380405
hamilt::folding_HR(*this->hR_tmp, this->hK->data(), this->kvec_d[ik], ncol, 0);
381406
}
382407

@@ -388,4 +413,4 @@ template class TDEkinetic<hamilt::OperatorLCAO<double, double>>;
388413
template class TDEkinetic<hamilt::OperatorLCAO<std::complex<double>, double>>;
389414
template class TDEkinetic<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>;
390415

391-
}
416+
}

‎source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
99
#include "module_basis/module_ao/ORB_table_phi.h"
1010
#include "module_basis/module_ao/ORB_gaunt_table.h"
11+
#include "module_hamilt_lcao/module_tddft/td_velocity.h"
1112

1213

1314
namespace hamilt
@@ -75,6 +76,8 @@ class TDEkinetic<OperatorLCAO<TK,TR>> : public OperatorLCAO<TK, TR>
7576
void calculate_HR(void);
7677
virtual void set_HR_fixed(void*)override;
7778

79+
TD_Velocity td_velocity;
80+
7881
private:
7982
const UnitCell* ucell = nullptr;
8083

@@ -109,6 +112,8 @@ class TDEkinetic<OperatorLCAO<TK,TR>> : public OperatorLCAO<TK, TR>
109112

110113
bool hR_tmp_done = false;
111114
bool allocated = false;
115+
bool output_hR_done = false;
116+
bool out_mat_R = false;
112117
};
113118

114119
} // namespace hamilt

‎source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -47,9 +47,9 @@ template <typename TK, typename TR>
4747
void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::init_td(void)
4848
{
4949
//calculate At in cartesian coorinates.
50-
double l_norm[3]={GlobalC::ucell.a1.norm() ,GlobalC::ucell.a2.norm() ,GlobalC::ucell.a3.norm()};
50+
double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()};
5151
double (&A)[3] = elecstate::H_TDDFT_pw::At;
52-
cart_At = -(GlobalC::ucell.a1*A[0]/l_norm[0] + GlobalC::ucell.a2*A[1]/l_norm[1] + GlobalC::ucell.a3*A[2]/l_norm[2]);
52+
cart_At = -(this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2]);
5353
std::cout << "cart_At: " << cart_At[0] << " " <<cart_At[1]<< " " << cart_At[2] << std::endl;
5454
}
5555
// initialize_HR()
@@ -372,7 +372,7 @@ void hamilt::TDNonlocal<hamilt::OperatorLCAO<TK, TR>>::contributeHk(int ik)
372372
template<>
373373
void hamilt::TDNonlocal<hamilt::OperatorLCAO<std::complex<double>, double>>::contributeHk(int ik)
374374
{
375-
if (GlobalV::ESOLVER_TYPE != "tddft" || elecstate::H_TDDFT_pw::stype != 1)
375+
if (TD_Velocity::tddft_velocity == false)
376376
{
377377
return;
378378
}

‎source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h"
99
#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
1010
#include "module_elecstate/potentials/H_TDDFT_pw.h"
11+
#include "module_hamilt_lcao/module_tddft/td_velocity.h"
1112

1213
namespace hamilt
1314
{

‎source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp

Lines changed: 79 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "spar_u.h"
44
#include "spar_exx.h"
55
#include "module_hamilt_lcao/module_hcontainer/hcontainer.h"
6+
#include "module_hamilt_lcao/module_tddft/td_velocity.h"
67

78
void sparse_format::cal_HSR(
89
const Parallel_Orbitals &pv,
@@ -25,12 +26,24 @@ void sparse_format::cal_HSR(
2526
hamilt::HamiltLCAO<std::complex<double>, double>* p_ham_lcao =
2627
dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(p_ham);
2728

28-
sparse_format::cal_HContainer_d(
29+
if(TD_Velocity::tddft_velocity)
30+
{
31+
sparse_format::cal_HContainer_td(
32+
pv,
33+
current_spin,
34+
sparse_thr,
35+
*(p_ham_lcao->getHR()),
36+
TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin]);
37+
}
38+
else
39+
{
40+
sparse_format::cal_HContainer_d(
2941
pv,
3042
current_spin,
3143
sparse_thr,
3244
*(p_ham_lcao->getHR()),
3345
lm.HR_sparse[current_spin]);
46+
}
3447

3548
sparse_format::cal_HContainer_d(
3649
pv,
@@ -202,6 +215,49 @@ void sparse_format::cal_HContainer_cd(
202215
return;
203216
}
204217

218+
void sparse_format::cal_HContainer_td(
219+
const Parallel_Orbitals &pv,
220+
const int &current_spin,
221+
const double &sparse_thr,
222+
const hamilt::HContainer<double>& hR,
223+
std::map<Abfs::Vector3_Order<int>,
224+
std::map<size_t, std::map<size_t, std::complex<double>>>>& target)
225+
{
226+
ModuleBase::TITLE("sparse_format","cal_HContainer_td");
227+
228+
auto row_indexes = pv.get_indexes_row();
229+
auto col_indexes = pv.get_indexes_col();
230+
for(int iap=0;iap<hR.size_atom_pairs();++iap)
231+
{
232+
int atom_i = hR.get_atom_pair(iap).get_atom_i();
233+
int atom_j = hR.get_atom_pair(iap).get_atom_j();
234+
int start_i = pv.atom_begin_row[atom_i];
235+
int start_j = pv.atom_begin_col[atom_j];
236+
int row_size = pv.get_row_size(atom_i);
237+
int col_size = pv.get_col_size(atom_j);
238+
for(int iR=0;iR<hR.get_atom_pair(iap).get_R_size();++iR)
239+
{
240+
auto& matrix = hR.get_atom_pair(iap).get_HR_values(iR);
241+
int* r_index = hR.get_atom_pair(iap).get_R_index(iR);
242+
Abfs::Vector3_Order<int> dR(r_index[0], r_index[1], r_index[2]);
243+
for(int i=0;i<row_size;++i)
244+
{
245+
int mu = row_indexes[start_i+i];
246+
for(int j=0;j<col_size;++j)
247+
{
248+
int nu = col_indexes[start_j+j];
249+
const auto& value_tmp = std::complex<double>(matrix.get_value(i,j) , 0.0);
250+
if(std::abs(value_tmp)>sparse_thr)
251+
{
252+
target[dR][mu][nu] += value_tmp;
253+
}
254+
}
255+
}
256+
}
257+
}
258+
259+
return;
260+
}
205261

206262
// in case there are elements smaller than the threshold
207263
void sparse_format::clear_zero_elements(
@@ -232,6 +288,28 @@ void sparse_format::clear_zero_elements(
232288
}
233289
}
234290
}
291+
if(TD_Velocity::tddft_velocity)
292+
{
293+
for (auto &R_loop : TD_Velocity::td_vel_op->HR_sparse_td_vel[current_spin])
294+
{
295+
for (auto &row_loop : R_loop.second)
296+
{
297+
auto &col_map = row_loop.second;
298+
auto iter = col_map.begin();
299+
while (iter != col_map.end())
300+
{
301+
if (std::abs(iter->second) <= sparse_thr)
302+
{
303+
col_map.erase(iter++);
304+
}
305+
else
306+
{
307+
iter++;
308+
}
309+
}
310+
}
311+
}
312+
}
235313

236314
for (auto &R_loop : lm.SR_sparse)
237315
{

‎source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,14 @@ namespace sparse_format
3131
std::map<Abfs::Vector3_Order<int>,
3232
std::map<size_t, std::map<size_t, std::complex<double>>>>& target);
3333

34+
void cal_HContainer_td(
35+
const Parallel_Orbitals &pv,
36+
const int &current_spin,
37+
const double &sparse_threshold,
38+
const hamilt::HContainer<double>& hR,
39+
std::map<Abfs::Vector3_Order<int>,
40+
std::map<size_t, std::map<size_t, std::complex<double>>>>& target);
41+
3442
void clear_zero_elements(
3543
LCAO_Matrix &lm,
3644
const int &current_spin,

‎source/module_hamilt_lcao/module_tddft/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ if(ENABLE_LCAO)
77
norm_psi.cpp
88
propagator.cpp
99
upsi.cpp
10+
td_velocity.cpp
1011
)
1112

1213
add_library(

0 commit comments

Comments
 (0)