Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions source/source_esolver/esolver_ks_lcao_tddft.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
// 1) before_scf (electronic iteration loops)
//----------------------------------------------------------------
this->before_scf(ucell, istep); // From ESolver_KS_LCAO

td_p->initialize_phase_hybrid(ucell, dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, TR>*>(this->p_hamilt)->getHR());
// Initialize the moving spatial gauge
if (use_td_moving_gauge && this->td_mg_ == nullptr)
{
Expand All @@ -113,7 +113,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)

if (PARAM.inp.td_stype == 2)
{
this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(),TD_info::cart_At);
}
else
{
Expand Down Expand Up @@ -594,7 +594,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::weight_dm_rho(const UnitCell& ucell)
elecstate::cal_dm_psi(this->dmat.dm->get_paraV_pointer(), this->pelec->wg, this->psi[0], *this->dmat.dm);
if (PARAM.inp.td_stype == 2)
{
this->dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
this->dmat.dm->cal_DMR_td(td_p->get_phase_hybrid(), TD_info::cart_At);
}
else
{
Expand Down
15 changes: 12 additions & 3 deletions source/source_estate/module_dm/cal_dm_psi.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -70,11 +70,11 @@ void cal_dm_psi(const Parallel_Orbitals* ParaV,

return;
}

template <typename TR>
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)
elecstate::DensityMatrix<std::complex<double>, TR>& DM)
{
ModuleBase::TITLE("elecstate", "cal_dm_psi");
ModuleBase::timer::start("elecstate", "cal_dm_psi");
Expand Down Expand Up @@ -268,5 +268,14 @@ void psiMulPsi(const psi::Psi<std::complex<double>>& psi1,
dm_out,
nlocal);
}

template
void cal_dm_psi(const Parallel_Orbitals* ParaV,
const ModuleBase::matrix& wg,
const psi::Psi<std::complex<double>>& wfc,
elecstate::DensityMatrix<std::complex<double>, std::complex<double>>& DM);
template
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);
} // namespace elecstate
3 changes: 2 additions & 1 deletion source/source_estate/module_dm/cal_dm_psi.h
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@ namespace elecstate
void cal_dm_psi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& wg, const psi::Psi<double>& wfc, elecstate::DensityMatrix<double, double>& DM);

// for Multi-k case where DMK is std::complex<double>
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);
template <typename TR>
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);

// for Gamma-Only case with MPI
void psiMulPsiMpi(const psi::Psi<double>& psi1,
Expand Down
22 changes: 12 additions & 10 deletions source/source_estate/module_dm/density_matrix.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -214,7 +214,7 @@ template <typename TK, typename TR_in, typename TR_out>
void DensityMatrix_Tools::cal_DMR_td(
const DensityMatrix<TK, TR_in> &dm,
std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
const UnitCell& ucell,
const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid,
const ModuleBase::Vector3<double> At,
const int ik_in)
{
Expand Down Expand Up @@ -262,19 +262,21 @@ void DensityMatrix_Tools::cal_DMR_td(
}
#endif
target_DMR_mat_vec[iR] = target_mat->get_pointer();
//cal tddft phase for hybrid gauge
const ModuleBase::Vector3<double> dtau = ucell.cal_dtau(iat1, iat2, R_index);
const double arg_td = At * dtau * ucell.lat0;
for(int ik = 0; ik < dm._nk; ++ik)
{
if(ik_in >= 0 && ik_in != ik) { continue; }
// cal k_phase
// if TK==std::complex<double>, kphase is e^{ikR}
const ModuleBase::Vector3<double> dR(R_index[0], R_index[1], R_index[2]);
const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI + arg_td;
const double arg = (dm._kvec_d[ik] * dR) * ModuleBase::TWO_PI;
double sinp, cosp;
ModuleBase::libm::sincos(arg, &sinp, &cosp);
kphase_vec[ik][iR] = TK(cosp, sinp);
if(PARAM.inp.td_stype==2)
{
//phase for hybrid gauge tddft
kphase_vec[ik][iR] *= phase_hybrid.at(R_index);
}
}
}

Expand Down Expand Up @@ -353,20 +355,20 @@ void DensityMatrix_Tools::cal_DMR_td(
ModuleBase::timer::end("DensityMatrix", "cal_DMR_td");
}
template <>
void DensityMatrix<double, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
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)
{
return;
}
template <>
void DensityMatrix<std::complex<double>, double>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
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)
{
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in);
}

template <>
void DensityMatrix<std::complex<double>, std::complex<double>>::cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in)
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)
{
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, ucell, At, ik_in);
DensityMatrix_Tools::cal_DMR_td(*this, this->_DMR, phase_hybrid, At, ik_in);
}


Expand Down
6 changes: 3 additions & 3 deletions source/source_estate/module_dm/density_matrix.h
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ namespace DensityMatrix_Tools
extern void cal_DMR_td(
const DensityMatrix<TK, TR_in> &dm,
std::vector<hamilt::HContainer<TR_out>*> &dmR_out,
const UnitCell& ucell,
const std::map<ModuleBase::Vector3<int>, std::complex<double>>& phase_hybrid,
const ModuleBase::Vector3<double> At,
const int ik_in);

Expand Down Expand Up @@ -225,7 +225,7 @@ class DensityMatrix
* if ik_in < 0, calculate all k-points
* if ik_in >= 0, calculate only one k-point without summing over k-points
*/
void cal_DMR_td(const UnitCell& ucell, const ModuleBase::Vector3<double> At, const int ik_in = -1);
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);

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

friend void DensityMatrix_Tools::cal_DMR<TK,TR>(const DensityMatrix<TK, TR> &dm, std::vector<hamilt::HContainer<TR>*> &dmR_out, const int ik_in);
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);
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);
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);
};

Expand Down
2 changes: 1 addition & 1 deletion source/source_estate/module_dm/init_dm.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ void elecstate::init_dm(UnitCell& ucell,
elecstate::cal_dm_psi(dmat.dm->get_paraV_pointer(), pelec->wg, *psi, *dmat.dm);
if (PARAM.inp.esolver_type!="tddft" && PARAM.inp.td_stype == 2)
{
dmat.dm->cal_DMR_td(ucell, TD_info::cart_At);
dmat.dm->cal_DMR_td(TD_info::td_vel_op->get_phase_hybrid(), TD_info::cart_At);
}
else
{
Expand Down
6 changes: 3 additions & 3 deletions source/source_io/module_ctrl/ctrl_output_td.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,16 +42,16 @@ void ctrl_output_td(const UnitCell& ucell,
{
if (TD_info::out_current_k)
{
ModuleIO::write_current_eachk<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
ModuleIO::write_current_eachk<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA);
}
else
{
ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, td_p, RA);
}
}
else if(TD_info::out_current==2)
{
ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
}
// (3) Output file for restart
if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0
Expand Down
Loading
Loading