Skip to content

Commit b16fcff

Browse files
ESROAMERmohanchen
ESROAMER
and
mohanchen
authored
Support using hybrid functionals for RT-TDDFT calculations. (#6886)
* Add files via upload * Add files via upload * Add files via upload * Add files via upload * Delete source/ctrl_output_td.h * Add files via upload * Add files via upload * Add files via upload * Add files via upload * Update td_info.cpp * Update td_current_io_comm.cpp --------- Co-authored-by: Mohan Chen <mohanchen@pku.edu.cn>
1 parent 272d6b4 commit b16fcff

30 files changed

Lines changed: 1395 additions & 97 deletions

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

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4133,11 +4133,12 @@ These variables are used to control berry phase and wannier90 interface paramete
41334133

41344134
### out_current
41354135

4136-
- **Type**: Boolean
4136+
- **Type**: Integer
41374137
- **Description**:
4138-
- True: Output current.
4139-
- False: Do not output current.
4140-
- **Default**: False
4138+
- 0: Do not output current.
4139+
- 1: Output current using the two-center integral, faster.
4140+
- 2: Output current using the matrix commutation, more precise.
4141+
- **Default**: 0
41414142

41424143
### out_current_k
41434144

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -553,6 +553,7 @@ OBJS_IO=input_conv.o\
553553
write_dipole.o\
554554
write_init.o\
555555
td_current_io.o\
556+
td_current_io_comm.o\
556557
write_libxc_r.o\
557558
output_log.o\
558559
output_mat_sparse.o\

‎source/source_esolver/esolver_ks_lcao_tddft.cpp‎

Lines changed: 15 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::before_all_runners(UnitCell& ucell, cons
5555
// Run before_all_runners in ESolver_KS_LCAO
5656
ESolver_KS_LCAO<std::complex<double>, TR>::before_all_runners(ucell, inp);
5757

58-
td_p = new TD_info(&ucell);
58+
td_p = new TD_info(&ucell, this->pv, this->orb_);
5959
TD_info::td_vel_op = td_p;
6060
totstep += TD_info::estep_shift;
6161

@@ -90,7 +90,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
9090
ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
9191

9292
// Initialize velocity operator for current calculation
93-
if (PARAM.inp.td_stype != 1 && TD_info::out_current)
93+
if (PARAM.inp.td_stype != 1 && TD_info::out_current == 1)
9494
{
9595
// initialize the velocity operator
9696
velocity_mat = new Velocity_op<TR>(&ucell,
@@ -203,7 +203,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
203203
}
204204
}
205205

206-
if (PARAM.inp.td_stype != 1 && TD_info::out_current)
206+
if(PARAM.inp.td_stype != 1 && TD_info::out_current == 1)
207207
{
208208
delete velocity_mat;
209209
}
@@ -296,6 +296,12 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
296296
srho.begin(is, this->chr, this->pw_rho, ucell.symm);
297297
}
298298
}
299+
#ifdef __EXX
300+
if (GlobalC::exx_info.info_ri.real_number)
301+
this->exx_nao.exd->exx_hamilt2rho(*this->pelec, this->pv, iter);
302+
else
303+
this->exx_nao.exc->exx_hamilt2rho(*this->pelec, this->pv, iter);
304+
#endif
299305

300306
// Calculate delta energy
301307
this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
@@ -474,6 +480,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::after_scf(UnitCell& ucell, const int ist
474480
std::cout << " Potential (Ry): " << std::setprecision(15) << this->pelec->f_en.etot << std::endl;
475481

476482
// Output dipole, current, etc.
483+
auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, TR>*>(this->p_hamilt);
477484
ModuleIO::ctrl_output_td<TR>(ucell,
478485
this->chr.rho_save,
479486
this->chr.rhopw,
@@ -485,8 +492,12 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::after_scf(UnitCell& ucell, const int ist
485492
&this->pv,
486493
this->orb_,
487494
this->velocity_mat,
495+
this->gd,
496+
hamilt_lcao,
488497
this->RA,
489-
this->td_p);
498+
this->td_p,
499+
this->exx_nao
500+
);
490501

491502
ModuleBase::timer::tick(this->classname, "after_scf");
492503
}

‎source/source_io/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,7 @@ list(APPEND objects
3434
write_init.cpp
3535
write_mlkedf_descriptors.cpp
3636
td_current_io.cpp
37+
td_current_io_comm.cpp
3738
write_libxc_r.cpp
3839
output_log.cpp
3940
para_json.cpp

‎source/source_io/ctrl_output_td.cpp‎

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -20,8 +20,12 @@ void ctrl_output_td(const UnitCell& ucell,
2020
const Parallel_Orbitals* pv,
2121
const LCAO_Orbitals& orb,
2222
const Velocity_op<TR>* velocity_mat,
23+
const Grid_Driver& grid,
24+
hamilt::HamiltLCAO<std::complex<double>, TR>* p_hamilt,
2325
Record_adj& RA,
24-
TD_info* td_p)
26+
TD_info* td_p,
27+
const Exx_NAO<std::complex<double>>& exx_nao
28+
)
2529
{
2630
ModuleBase::TITLE("ModuleIO", "ctrl_output_td");
2731

@@ -46,7 +50,7 @@ void ctrl_output_td(const UnitCell& ucell,
4650
ModuleBase::WARNING_QUIT("ModuleIO::ctrl_output_td", "Failed to cast ElecState to ElecStateLCAO");
4751
}
4852

49-
if (TD_info::out_current)
53+
if (TD_info::out_current == 1)
5054
{
5155
if (TD_info::out_current_k)
5256
{
@@ -57,6 +61,10 @@ void ctrl_output_td(const UnitCell& ucell,
5761
ModuleIO::write_current<TR>(ucell, istep, psi, pelec, kv, intor, pv, orb, velocity_mat, RA);
5862
}
5963
}
64+
else if(TD_info::out_current==2)
65+
{
66+
ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
67+
}
6068

6169
// (3) Output file for restart
6270
if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0
@@ -88,8 +96,12 @@ template void ctrl_output_td<double>(const UnitCell&,
8896
const Parallel_Orbitals*,
8997
const LCAO_Orbitals&,
9098
const Velocity_op<double>*,
99+
const Grid_Driver&,
100+
hamilt::HamiltLCAO<std::complex<double>, double>*,
91101
Record_adj&,
92-
TD_info*);
102+
TD_info*,
103+
const Exx_NAO<std::complex<double>>&
104+
);
93105

94106
template void ctrl_output_td<std::complex<double>>(const UnitCell&,
95107
double**,
@@ -102,7 +114,11 @@ template void ctrl_output_td<std::complex<double>>(const UnitCell&,
102114
const Parallel_Orbitals*,
103115
const LCAO_Orbitals&,
104116
const Velocity_op<std::complex<double>>*,
117+
const Grid_Driver&,
118+
hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*,
105119
Record_adj&,
106-
TD_info*);
120+
TD_info*,
121+
const Exx_NAO<std::complex<double>>&
122+
);
107123

108124
} // namespace ModuleIO

‎source/source_io/ctrl_output_td.h‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,11 @@
1111
#include "source_lcao/module_rt/velocity_op.h"
1212
#include "source_lcao/record_adj.h"
1313
#include "source_psi/psi.h"
14+
#include "source_lcao/hamilt_lcao.h"
15+
#include "source_lcao/setup_exx.h"
16+
#ifdef __EXX
17+
#include <RI/global/Tensor.h>
18+
#endif
1419

1520
namespace ModuleIO
1621
{
@@ -27,8 +32,12 @@ void ctrl_output_td(const UnitCell& ucell,
2732
const Parallel_Orbitals* pv,
2833
const LCAO_Orbitals& orb,
2934
const Velocity_op<TR>* velocity_mat,
35+
const Grid_Driver& grid,
36+
hamilt::HamiltLCAO<std::complex<double>, TR>* p_hamilt,
3037
Record_adj& RA,
31-
TD_info* td_p);
38+
TD_info* td_p,
39+
const Exx_NAO<std::complex<double>>& exx_nao
40+
);
3241

3342
} // namespace ModuleIO
3443

‎source/source_io/module_parameter/input_parameter.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -403,7 +403,7 @@ struct Input_para
403403
int out_wfc_lcao = 0; ///< output the wave functions in local basis.
404404
bool out_dipole = false; ///< output the dipole or not
405405
bool out_efield = false; ///< output the efield or not
406-
bool out_current = false; ///< output the current or not
406+
int out_current = 0; ///< output the current or not
407407
bool out_current_k = false; ///< output tddft current for all k points
408408
bool out_vecpot = false; ///< output the vector potential or not
409409
bool restart_save = false; ///< restart //Peize Lin add 2020-04-04

‎source/source_io/read_input_item_exx_dftu.cpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,13 @@ void ReadInput::item_exx()
164164
Input_Item item("exx_separate_loop");
165165
item.annotation = "if 1, a two-step method is employed, else it will "
166166
"start with a GGA-Loop, and then Hybrid-Loop";
167+
item.reset_value = [](const Input_Item& item, Parameter& para) {
168+
if (para.input.esolver_type == "tddft" && para.input.exx_separate_loop)
169+
{
170+
GlobalV::ofs_running << "For RT-TDDFT with hybrid functionals, only exx_separate_loop = 0 is supported" << std::endl;
171+
para.input.exx_separate_loop = false;
172+
}
173+
};
167174
read_sync_bool(input.exx_separate_loop);
168175
this->add_item(item);
169176
}

‎source/source_io/read_input_item_output.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -443,7 +443,7 @@ void ReadInput::item_output()
443443
{
444444
Input_Item item("out_current");
445445
item.annotation = "output current or not";
446-
read_sync_bool(input.out_current);
446+
read_sync_int(input.out_current);
447447
this->add_item(item);
448448
}
449449
{

‎source/source_io/td_current_io.h‎

Lines changed: 74 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,10 @@
66
#include "source_estate/module_dm/density_matrix.h"
77
#include "source_psi/psi.h"
88
#include "source_lcao/module_rt/velocity_op.h"
9+
#include "source_lcao/setup_exx.h"
10+
#ifdef __EXX
11+
#include <RI/global/Tensor.h>
12+
#endif
913

1014
namespace ModuleIO
1115
{
@@ -33,7 +37,22 @@ void write_current(const UnitCell& ucell,
3337
const LCAO_Orbitals& orb,
3438
const Velocity_op<TR>* cal_current,
3539
Record_adj& ra);
36-
40+
/// @brief func to output current calculated using i[r,H] directly
41+
template <typename TR>
42+
void write_current(
43+
const UnitCell& ucell,
44+
const Grid_Driver& GridD,
45+
const int istep,
46+
const psi::Psi<std::complex<double>>* psi,
47+
const elecstate::ElecState* pelec,
48+
const K_Vectors& kv,
49+
const Parallel_Orbitals* pv,
50+
const LCAO_Orbitals& orb,
51+
cal_r_overlap_R& r_calculator,
52+
const hamilt::HContainer<TR>* sR,
53+
const hamilt::HContainer<TR>* hR,
54+
const Exx_NAO<std::complex<double>>& exx_nao
55+
);
3756
/// @brief calculate sum_n[𝜌_(𝑛𝑘,𝜇𝜈)] for current calculation
3857
void cal_tmp_DM_k(const UnitCell& ucell,
3958
elecstate::DensityMatrix<std::complex<double>, double>& DM_real,
@@ -47,7 +66,61 @@ void cal_tmp_DM(const UnitCell& ucell,
4766
elecstate::DensityMatrix<std::complex<double>, double>& DM_real,
4867
elecstate::DensityMatrix<std::complex<double>, double>& DM_imag,
4968
const int nspin);
69+
void set_rR_from_hR(const UnitCell& ucell,
70+
const Grid_Driver& GridD,
71+
const LCAO_Orbitals& orb,
72+
const Parallel_Orbitals* pv,
73+
cal_r_overlap_R& r_calculator,
74+
const hamilt::HContainer<std::complex<double>>* hR,
75+
ModuleBase::Vector3<hamilt::HContainer<double>*>& rR);
76+
template <typename TR>
77+
void sum_HR(
78+
const UnitCell& ucell,
79+
const Parallel_Orbitals& pv,
80+
const K_Vectors& kv,
81+
const hamilt::HContainer<TR>* hR,
82+
hamilt::HContainer<std::complex<double>>* full_hR,
83+
const Exx_NAO<std::complex<double>>& exx_nao
84+
);
85+
86+
template <typename Tadd, typename Tfull>
87+
void add_HR(const hamilt::HContainer<Tadd>* hR, hamilt::HContainer<Tfull>* full_hR);
88+
89+
void init_from_adj(const UnitCell& ucell,
90+
const Grid_Driver& GridD,
91+
const LCAO_Orbitals& orb,
92+
const Parallel_Orbitals* pv,
93+
std::vector<AdjacentAtomInfo>& adjs_all,
94+
ModuleBase::Vector3<hamilt::HContainer<double>*>& rR);
95+
template <typename TR, typename TA>
96+
void init_from_hR(const hamilt::HContainer<TR>* hR, hamilt::HContainer<TA>* aimR);
97+
template <typename TR>
98+
void cal_velocity_basis_k(const UnitCell& ucell,
99+
const LCAO_Orbitals& orb,
100+
const Parallel_Orbitals* pv,
101+
const K_Vectors& kv,
102+
const ModuleBase::Vector3<hamilt::HContainer<double>*>& rR,
103+
const hamilt::HContainer<TR>& sR,
104+
const hamilt::HContainer<std::complex<double>>& hR,
105+
std::vector<ModuleBase::Vector3<std::complex<double>*>>& velocity_basis_k);
50106

107+
void cal_velocity_matrix(const psi::Psi<std::complex<double>>* psi,
108+
const Parallel_Orbitals* pv,
109+
const K_Vectors& kv,
110+
const std::vector<ModuleBase::Vector3<std::complex<double>*>>& velocity_basis_k,
111+
std::vector<std::array<ModuleBase::ComplexMatrix, 3>>& velocity_k);
112+
template <typename TR>
113+
void cal_current_comm_k(const UnitCell& ucell,
114+
const Grid_Driver& GridD,
115+
const LCAO_Orbitals& orb,
116+
const Parallel_Orbitals* pv,
117+
const K_Vectors& kv,
118+
cal_r_overlap_R& r_calculator,
119+
const hamilt::HContainer<TR>& sR,
120+
const hamilt::HContainer<std::complex<double>>& hR,
121+
const psi::Psi<std::complex<double>>* psi,
122+
const elecstate::ElecState* pelec,
123+
std::vector<ModuleBase::Vector3<double>>& current_k);
51124
#endif // __LCAO
52125
} // namespace ModuleIO
53126
#endif

0 commit comments

Comments
 (0)