Skip to content

Commit fa7d4a4

Browse files
AsTonyshmentAsTonyshment
authored
Fix wrong output of S(k) in TDDFT calculation (#4148)
* Fix wrong Sk ouput in rt-TDDFT * Did nothing, just formatting with clang-format * Add const qualifiers to parameters in hamilt2density to prevent modifications
1 parent a2798b2 commit fa7d4a4

File tree

1 file changed

+83
-102
lines changed

1 file changed

+83
-102
lines changed

‎source/module_esolver/esolver_ks_lcao_tddft.cpp

Lines changed: 83 additions & 102 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,9 @@
44
#include "module_io/dipole_io.h"
55
#include "module_io/dm_io.h"
66
#include "module_io/rho_io.h"
7+
#include "module_io/td_current_io.h"
78
#include "module_io/write_HS.h"
89
#include "module_io/write_HS_R.h"
9-
#include "module_io/td_current_io.h"
1010

1111
//--------------temporary----------------------------
1212
#include "module_base/blas_connector.h"
@@ -39,7 +39,6 @@ ESolver_KS_LCAO_TDDFT::ESolver_KS_LCAO_TDDFT()
3939
basisname = "LCAO";
4040
}
4141

42-
4342
ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT()
4443
{
4544
// this->orb_con.clear_after_ions(GlobalC::UOT, GlobalC::ORB, GlobalV::deepks_setorb, GlobalC::ucell.infoNL.nproj);
@@ -79,8 +78,8 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
7978
this->pelec = new elecstate::ElecStateLCAO_TDDFT(&(this->chr),
8079
&(kv),
8180
kv.nks,
82-
&(this->LOC),
83-
&(this->GK), // mohan add 2024-04-01
81+
&(this->LOC),
82+
&(this->GK), // mohan add 2024-04-01
8483
&(this->LOWF),
8584
this->pw_rho,
8685
pw_big);
@@ -105,7 +104,8 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
105104
this->LOC.ParaV = this->LOWF.ParaV = this->LM.ParaV;
106105

107106
// init DensityMatrix
108-
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->init_DM(&kv, this->LM.ParaV, GlobalV::NSPIN);
107+
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)
108+
->init_DM(&kv, this->LM.ParaV, GlobalV::NSPIN);
109109

110110
// init Psi, HSolver, ElecState, Hamilt
111111
if (this->phsol == nullptr)
@@ -130,12 +130,8 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
130130
this->pelec_td = dynamic_cast<elecstate::ElecStateLCAO_TDDFT*>(this->pelec);
131131
}
132132

133-
void ESolver_KS_LCAO_TDDFT::hamilt2density(
134-
int istep,
135-
int iter,
136-
double ethr)
133+
void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr)
137134
{
138-
139135
pelec->charge->save_rho_before_sum_band();
140136

141137
if (wf.init_wfc == "file")
@@ -184,12 +180,8 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(
184180
this->pelec->f_en.demet = 0.0;
185181
if (this->psi != nullptr)
186182
{
187-
this->phsol->solve(
188-
this->p_hamilt,
189-
this->psi[0],
190-
this->pelec_td,
191-
GlobalV::KS_SOLVER);
192-
}
183+
this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec_td, GlobalV::KS_SOLVER);
184+
}
193185
}
194186
else
195187
{
@@ -211,13 +203,9 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(
211203
for (int ib = 0; ib < GlobalV::NBANDS; ib++)
212204
{
213205
std::setprecision(6);
214-
GlobalV::ofs_running << ik + 1
215-
<< " "
216-
<< ib + 1
217-
<< " "
218-
<< this->pelec_td->wg(ik, ib)
219-
<< std::endl;
220-
}
206+
GlobalV::ofs_running << ik + 1 << " " << ib + 1 << " " << this->pelec_td->wg(ik, ib)
207+
<< std::endl;
208+
}
221209
}
222210
GlobalV::ofs_running << std::endl;
223211
GlobalV::ofs_running
@@ -239,12 +227,8 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(
239227
Symmetry_rho srho;
240228
for (int is = 0; is < GlobalV::NSPIN; is++)
241229
{
242-
srho.begin(is,
243-
*(pelec->charge),
244-
pw_rho,
245-
GlobalC::Pgrid,
246-
GlobalC::ucell.symm);
247-
}
230+
srho.begin(is, *(pelec->charge), pw_rho, GlobalC::Pgrid, GlobalC::ucell.symm);
231+
}
248232
}
249233

250234
// (6) compute magnetization, only for spin==2
@@ -280,29 +264,29 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
280264
this->p_hamilt->matrix(h_mat, s_mat);
281265
if (hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs[0])
282266
{
283-
ModuleIO::save_mat(istep,
284-
h_mat.p,
285-
GlobalV::NLOCAL,
286-
bit,
287-
hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs[1],
288-
1,
289-
GlobalV::out_app_flag,
290-
"H",
291-
"data-" + std::to_string(ik),
292-
*this->LOWF.ParaV,
293-
GlobalV::DRANK);
294-
295-
ModuleIO::save_mat(istep,
296-
h_mat.p,
297-
GlobalV::NLOCAL,
298-
bit,
299-
hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs[1],
300-
1,
301-
GlobalV::out_app_flag,
302-
"S",
303-
"data-" + std::to_string(ik),
304-
*this->LOWF.ParaV,
305-
GlobalV::DRANK);
267+
ModuleIO::save_mat(istep,
268+
h_mat.p,
269+
GlobalV::NLOCAL,
270+
bit,
271+
hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs[1],
272+
1,
273+
GlobalV::out_app_flag,
274+
"H",
275+
"data-" + std::to_string(ik),
276+
*this->LOWF.ParaV,
277+
GlobalV::DRANK);
278+
279+
ModuleIO::save_mat(istep,
280+
s_mat.p,
281+
GlobalV::NLOCAL,
282+
bit,
283+
hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs[1],
284+
1,
285+
GlobalV::out_app_flag,
286+
"S",
287+
"data-" + std::to_string(ik),
288+
*this->LOWF.ParaV,
289+
GlobalV::DRANK);
306290
}
307291
}
308292
}
@@ -313,14 +297,14 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
313297
if (elecstate::ElecStateLCAO<std::complex<double>>::out_wfc_lcao)
314298
{
315299
elecstate::ElecStateLCAO<std::complex<double>>::out_wfc_flag
316-
= elecstate::ElecStateLCAO<std::complex<double>>::out_wfc_lcao;
300+
= elecstate::ElecStateLCAO<std::complex<double>>::out_wfc_lcao;
317301
}
318302
for (int ik = 0; ik < kv.nks; ik++)
319303
{
320304
if (istep % GlobalV::out_interval == 0)
321305
{
322-
this->psi[0].fix_k(ik);
323-
this->pelec->print_psi(this->psi[0], istep);
306+
this->psi[0].fix_k(ik);
307+
this->pelec->print_psi(this->psi[0], istep);
324308
}
325309
}
326310
elecstate::ElecStateLCAO<std::complex<double>>::out_wfc_flag = 0;
@@ -329,11 +313,11 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
329313
// Calculate new potential according to new Charge Density
330314
if (!this->conv_elec)
331315
{
332-
if (GlobalV::NSPIN == 4)
333-
{
334-
GlobalC::ucell.cal_ux();
335-
}
336-
this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell);
316+
if (GlobalV::NSPIN == 4)
317+
{
318+
GlobalC::ucell.cal_ux();
319+
}
320+
this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell);
337321
this->pelec->f_en.descf = this->pelec->cal_delta_escf();
338322
}
339323
else
@@ -344,19 +328,19 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
344328
// store wfc and Hk laststep
345329
if (istep >= (wf.init_wfc == "file" ? 0 : 1) && this->conv_elec)
346330
{
347-
if (this->psi_laststep == nullptr)
348-
{
331+
if (this->psi_laststep == nullptr)
332+
{
349333
#ifdef __MPI
350-
this->psi_laststep = new psi::Psi<std::complex<double>>(kv.nks,
351-
this->LOWF.ParaV->ncol_bands,
352-
this->LOWF.ParaV->nrow,
353-
nullptr);
334+
this->psi_laststep = new psi::Psi<std::complex<double>>(kv.nks,
335+
this->LOWF.ParaV->ncol_bands,
336+
this->LOWF.ParaV->nrow,
337+
nullptr);
354338
#else
355-
this->psi_laststep = new psi::Psi<std::complex<double>>(kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr);
339+
this->psi_laststep = new psi::Psi<std::complex<double>>(kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr);
356340
#endif
357-
}
341+
}
358342

359-
if (td_htype == 1)
343+
if (td_htype == 1)
360344
{
361345
if (this->Hk_laststep == nullptr)
362346
{
@@ -383,10 +367,10 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
383367
this->psi->fix_k(ik);
384368
this->psi_laststep->fix_k(ik);
385369
int size0 = psi->get_nbands() * psi->get_nbasis();
386-
for (int index = 0; index < size0; ++index)
387-
{
388-
psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index];
389-
}
370+
for (int index = 0; index < size0; ++index)
371+
{
372+
psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index];
373+
}
390374

391375
// store Hamiltonian
392376
if (td_htype == 1)
@@ -400,9 +384,8 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
400384
}
401385

402386
// calculate energy density matrix for tddft
403-
if (istep >= (wf.init_wfc == "file" ? 0 : 2)
404-
&& module_tddft::Evolve_elec::td_edm == 0)
405-
{
387+
if (istep >= (wf.init_wfc == "file" ? 0 : 2) && module_tddft::Evolve_elec::td_edm == 0)
388+
{
406389
this->cal_edm_tddft();
407390
}
408391
}
@@ -433,7 +416,6 @@ void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)
433416
}
434417
}
435418

436-
437419
void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
438420
{
439421
for (int is = 0; is < GlobalV::NSPIN; is++)
@@ -445,42 +427,41 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep)
445427
ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str());
446428
}
447429
}
448-
if(module_tddft::Evolve_elec::out_current == 1)
430+
if (module_tddft::Evolve_elec::out_current == 1)
449431
{
450-
elecstate::DensityMatrix<std::complex<double>, double>* tmp_DM =
451-
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM();
432+
elecstate::DensityMatrix<std::complex<double>, double>* tmp_DM
433+
= dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM();
452434
ModuleIO::write_current(istep,
453-
this->psi,
454-
pelec,
455-
kv,
456-
tmp_DM->get_paraV_pointer(),
457-
this->RA,
458-
this->LM, // mohan add 2024-04-02
459-
this->gen_h); // mohan add 2024-02
460-
}
435+
this->psi,
436+
pelec,
437+
kv,
438+
tmp_DM->get_paraV_pointer(),
439+
this->RA,
440+
this->LM, // mohan add 2024-04-02
441+
this->gen_h); // mohan add 2024-02
442+
}
461443
ESolver_KS_LCAO<std::complex<double>, double>::after_scf(istep);
462444
}
463445

464-
465446
// use the original formula (Hamiltonian matrix) to calculate energy density matrix
466447
void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void)
467448
{
468449
// mohan add 2024-03-27
469-
const int nlocal = GlobalV::NLOCAL;
470-
assert(nlocal>=0);
450+
const int nlocal = GlobalV::NLOCAL;
451+
assert(nlocal >= 0);
471452

472-
//this->LOC.edm_k_tddft.resize(kv.nks);
453+
// this->LOC.edm_k_tddft.resize(kv.nks);
473454
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->EDMK.resize(kv.nks);
474455
for (int ik = 0; ik < kv.nks; ++ik)
475456
{
476-
std::complex<double>* tmp_dmk =
477-
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->get_DMK_pointer(ik);
457+
std::complex<double>* tmp_dmk
458+
= dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->get_DMK_pointer(ik);
478459

479-
ModuleBase::ComplexMatrix& tmp_edmk =
480-
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->EDMK[ik];
460+
ModuleBase::ComplexMatrix& tmp_edmk
461+
= dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->EDMK[ik];
481462

482-
const Parallel_Orbitals* tmp_pv =
483-
dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->get_paraV_pointer();
463+
const Parallel_Orbitals* tmp_pv
464+
= dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM()->get_paraV_pointer();
484465

485466
#ifdef __MPI
486467

@@ -489,7 +470,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void)
489470
//! whether the long type is safe, needs more discussion
490471
const long nloc = this->LOC.ParaV->nloc;
491472

492-
//this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
473+
// this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
493474
tmp_edmk.create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
494475
complex<double>* Htmp = new complex<double>[nloc];
495476
complex<double>* Sinv = new complex<double>[nloc];
@@ -651,7 +632,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void)
651632
&one_int,
652633
this->LOC.ParaV->desc);
653634
zcopy_(&nloc, tmp4, &inc, tmp_edmk.c, &inc);
654-
//zcopy_(&nloc, tmp4, &inc, this->LOC.edm_k_tddft[ik].c, &inc);
635+
// zcopy_(&nloc, tmp4, &inc, this->LOC.edm_k_tddft[ik].c, &inc);
655636
delete[] Htmp;
656637
delete[] Sinv;
657638
delete[] tmp1;
@@ -661,7 +642,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void)
661642
delete[] ipiv;
662643
#else
663644
// for serial version
664-
//this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
645+
// this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
665646
tmp_edmk.create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow);
666647
ModuleBase::ComplexMatrix Sinv(nlocal, nlocal);
667648
ModuleBase::ComplexMatrix Htmp(nlocal, nlocal);
@@ -679,7 +660,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft(void)
679660
Sinv(i, j) = s_mat.p[i * nlocal + j];
680661
}
681662
}
682-
int INFO=0;
663+
int INFO = 0;
683664

684665
int lwork = 3 * nlocal - 1; // tmp
685666
std::complex<double>* work = new std::complex<double>[lwork];

0 commit comments

Comments
 (0)