Skip to content

Commit 260a811

Browse files
ESROAMERmohanchen
ESROAMER
and
mohanchen
authored
Upload simplified hybrid gauge rt-tddft (#7429)
* upload simplified hybrid gauge rt-tddft * Add files via upload * Initialize TD_info with phase_hybrid and static instance * Update current_tot.txt.ref * Update current_tot.txt.ref * Fix typo in README for test instructions * Update reference values in current_tot.txt.ref * Update result.ref --------- Co-authored-by: Mohan Chen <mohanchen@pku.edu.cn>
1 parent 4c1e6fa commit 260a811

27 files changed

Lines changed: 248 additions & 398 deletions

File tree

‎source/source_esolver/esolver_ks_lcao_tddft.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -101,7 +101,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
101101
// 1) before_scf (electronic iteration loops)
102102
//----------------------------------------------------------------
103103
this->before_scf(ucell, istep); // From ESolver_KS_LCAO
104-
104+
td_p->initialize_phase_hybrid(ucell, dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, TR>*>(this->p_hamilt)->getHR());
105105
// Initialize the moving spatial gauge
106106
if (use_td_moving_gauge && this->td_mg_ == nullptr)
107107
{
@@ -113,7 +113,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
113113

114114
if (PARAM.inp.td_stype == 2)
115115
{
116-
this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
116+
this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(),TD_info::cart_At);
117117
}
118118
else
119119
{
@@ -594,7 +594,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::weight_dm_rho(const UnitCell& ucell)
594594
elecstate::cal_dm_psi(this->dmat.dm->get_paraV_pointer(), this->pelec->wg, this->psi[0], *this->dmat.dm);
595595
if (PARAM.inp.td_stype == 2)
596596
{
597-
this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
597+
this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(), TD_info::cart_At);
598598
}
599599
else
600600
{

‎source/source_estate/module_dm/cal_dm_psi.cpp‎

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -70,11 +70,11 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,
7070

7171
return;
7272
}
73-
73+
template <typename TR>
7474
void cal_dm_psi(const Parallel_Orbitals* ParaV,
7575
const ModuleBase::matrix& wg,
7676
const psi::Psi<std::complex<double>>& wfc,
77-
elecstate::DensityMatrix<std::complex<double>, double>& DM)
77+
elecstate::DensityMatrix<std::complex<double>, TR>& DM)
7878
{
7979
ModuleBase::TITLE("elecstate", "cal_dm_psi");
8080
ModuleBase::timer::start("elecstate", "cal_dm_psi");
@@ -268,5 +268,14 @@ void psiMulPsi(const psi::Psi<std::complex<double>>& psi1,
268268
dm_out,
269269
nlocal);
270270
}
271-
271+
template
272+
void cal_dm_psi(const Parallel_Orbitals* ParaV,
273+
const ModuleBase::matrix& wg,
274+
const psi::Psi<std::complex<double>>& wfc,
275+
elecstate::DensityMatrix<std::complex<double>, std::complex<double>>& DM);
276+
template
277+
void cal_dm_psi(const Parallel_Orbitals* ParaV,
278+
const ModuleBase::matrix& wg,
279+
const psi::Psi<std::complex<double>>& wfc,
280+
elecstate::DensityMatrix<std::complex<double>, double>& DM);
272281
} // namespace elecstate

‎source/source_estate/module_dm/cal_dm_psi.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,8 @@ namespace elecstate
1111
void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<double>& wfc, elecstate::DensityMatrix<double, double>& DM);
1212

1313
// for Multi-k case where DMK is std::complex<double>
14-
void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<std::complex<double>>& wfc, elecstate::DensityMatrix<std::complex<double>, double>& DM);
14+
template <typename TR>
15+
void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<std::complex<double>>& wfc, elecstate::DensityMatrix<std::complex<double>, TR>& DM);
1516

1617
// for Gamma-Only case with MPI
1718
void psiMulPsiMpi(const psi::Psi<double>& psi1,

‎source/source_estate/module_dm/density_matrix.cpp‎

Lines changed: 12 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -214,7 +214,7 @@ template <typename TK, typename TR_in, typename TR_out>
214214
void DensityMatrix_Tools::cal_DMR_td(
215215
const DensityMatrix<TK, TR_in> &dm,
216216
std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
217-
const UnitCell& ucell,
217+
const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid,
218218
const ModuleBase::Vector3<double> At,
219219
const int ik_in)
220220
{
@@ -262,19 +262,21 @@ void DensityMatrix_Tools::cal_DMR_td(
262262
}
263263
#endif
264264
target_DMR_mat_vec[iR] = target_mat->get_pointer();
265-
//cal tddft phase for hybrid gauge
266-
const ModuleBase::Vector3<double> dtau = ucell.cal_dtau(iat1, iat2, R_index);
267-
const double arg_td = At * dtau * ucell.lat0;
268265
for(int ik = 0; ik < dm._nk; ++ik)
269266
{
270267
if(ik_in >= 0 && ik_in != ik) { continue; }
271268
// cal k_phase
272269
// if TK==std::complex<double>, kphase is e^{ikR}
273270
const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
274-
const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI + arg_td;
271+
const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
275272
double sinp, cosp;
276273
ModuleBase::libm::sincos(arg, &sinp, &cosp);
277274
kphase_vec[ik][iR] = TK(cosp, sinp);
275+
if(PARAM.inp.td_stype==2)
276+
{
277+
//phase for hybrid gauge tddft
278+
kphase_vec[ik][iR] *= phase_hybrid.at(R_index);
279+
}
278280
}
279281
}
280282

@@ -353,20 +355,20 @@ void DensityMatrix_Tools::cal_DMR_td(
353355
ModuleBase::timer::end("DensityMatrix", "cal_DMR_td");
354356
}
355357
template <>
356-
void DensityMatrix<double, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
358+
void DensityMatrix<double, double>::cal_DMR_td(const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid, const ModuleBase::Vector3<double> At, const int ik_in)
357359
{
358360
return;
359361
}
360362
template <>
361-
void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
363+
void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid, const ModuleBase::Vector3<double> At, const int ik_in)
362364
{
363-
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
365+
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in);
364366
}
365367

366368
template <>
367-
void DensityMatrix<std::complex<double>, std::complex<double>>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
369+
void DensityMatrix<std::complex<double>, std::complex<double>>::cal_DMR_td(const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid, const ModuleBase::Vector3<double> At, const int ik_in)
368370
{
369-
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
371+
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in);
370372
}
371373

372374

‎source/source_estate/module_dm/density_matrix.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ namespace DensityMatrix_Tools
4747
extern void cal_DMR_td(
4848
const DensityMatrix<TK, TR_in> &dm,
4949
std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
50-
const UnitCell& ucell,
50+
const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid,
5151
const ModuleBase::Vector3<double> At,
5252
const int ik_in);
5353

@@ -225,7 +225,7 @@ class DensityMatrix
225225
* if ik_in < 0, calculate all k-points
226226
* if ik_in >= 0, calculate only one k-point without summing over k-points
227227
*/
228-
void cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in = -1);
228+
void cal_DMR_td(const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid, const ModuleBase::Vector3<double> At, const int ik_in = -1);
229229

230230
/**
231231
* @brief calculate complex density matrix DMR with both real and imaginary part for noncollinear-spin calculation
@@ -327,7 +327,7 @@ class DensityMatrix
327327
TR* dmr_tmp_ = nullptr;
328328

329329
friend void DensityMatrix_Tools::cal_DMR<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const int ik_in);
330-
friend void DensityMatrix_Tools::cal_DMR_td<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in);
330+
friend void DensityMatrix_Tools::cal_DMR_td<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid, const ModuleBase::Vector3<double> At, const int ik_in);
331331
friend void DensityMatrix_Tools::cal_DMR_full<TK,TR>(const DensityMatrix<TK, TR> &dm, hamilt::HContainer<std::complex<double>>* dmR_out, const int ik_in);
332332
};
333333

‎source/source_estate/module_dm/init_dm.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ void elecstate::init_dm(UnitCell& ucell,
2525
elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm);
2626
if (PARAM.inp.esolver_type!="tddft" && PARAM.inp.td_stype == 2)
2727
{
28-
dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
28+
dmat.dm->cal_DMR_td(TD_info::td_vel_op->get_phase_hybrid(), TD_info::cart_At);
2929
}
3030
else
3131
{

‎source/source_io/module_ctrl/ctrl_output_td.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -42,16 +42,16 @@ void ctrl_output_td(const UnitCell& ucell,
4242
{
4343
if (TD_info::out_current_k)
4444
{
45-
ModuleIO::write_current_eachk<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
45+
ModuleIO::write_current_eachk<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA);
4646
}
4747
else
4848
{
49-
ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
49+
ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA);
5050
}
5151
}
5252
else if(TD_info::out_current==2)
5353
{
54-
ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
54+
ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
5555
}
5656
// (3) Output file for restart
5757
if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0

0 commit comments

Comments
 (0)