Skip to content

Commit

Permalink
fix bug in OperatorLREXX
Browse files Browse the repository at this point in the history
  • Loading branch information
@maki49
maki49 committed Aug 28, 2024
1 parent 8906098 commit 28fa245
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 12 deletions.
6 changes: 6 additions & 0 deletions source/module_lr/esolver_lrtd_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,9 @@ LR::ESolver_LR<T, TR>::ESolver_LR(ModuleESolver::ESolver_KS_LCAO<T, TR>&& ks_sol
this->move_exx_lri(ks_sol.exx_lri_complex);
} else // construct C, V from scratch
{
// set ccp_type according to the xc_kernel
if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; }
else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; }
this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
this->exx_lri->init(MPI_COMM_WORLD, this->kv);
this->exx_lri->cal_exx_ions(input.out_ri_cv);
Expand Down Expand Up @@ -362,6 +365,9 @@ LR::ESolver_LR<T, TR>::ESolver_LR(const Input_para& inp, UnitCell& ucell) : inpu
#ifdef __EXX
if ((xc_kernel == "hf" || xc_kernel == "hse") && this->input.lr_solver != "spectrum")
{
// set ccp_type according to the xc_kernel
if (xc_kernel == "hf") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; }
else if (xc_kernel == "hse") { exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; }
this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
this->exx_lri->init(MPI_COMM_WORLD, this->kv);
this->exx_lri->cal_exx_ions(input.out_ri_cv);
Expand Down
2 changes: 1 addition & 1 deletion source/module_lr/esolver_lrtd_lcao.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ namespace LR
std::shared_ptr<Exx_LRI<T>> exx_lri = nullptr;
void move_exx_lri(std::shared_ptr<Exx_LRI<double>>&);
void move_exx_lri(std::shared_ptr<Exx_LRI<std::complex<double>>>&);
const Exx_Info& exx_info;
Exx_Info& exx_info;
#endif
};
}
35 changes: 26 additions & 9 deletions source/module_lr/operator_casida/operator_lr_exx.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include "operator_lr_exx.h"
#include "module_lr/dm_trans/dm_trans.h"
#include "module_lr/utils/lr_util.h"
#include "module_lr/utils/lr_util_print.h"
namespace LR
{
template<typename T>
Expand Down Expand Up @@ -38,7 +39,7 @@ namespace LR
auto& D2d = this->Ds_onebase[is][iat1][std::make_pair(iat2, cell)];
for (int iw1 = 0;iw1 < ucell.atoms[it1].nw;++iw1)
for (int iw2 = 0;iw2 < ucell.atoms[it2].nw;++iw2)
D2d(iw1, iw2) = this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1)) * this->psi_ks_full(ik, iv, ucell.itiaiw2iwt(it2, ia2, iw2));
D2d(iw1, iw2) = this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1)) * this->psi_ks_full(ik, nocc + iv, ucell.itiaiw2iwt(it2, ia2, iw2));
}
}
}
Expand All @@ -61,7 +62,7 @@ namespace LR
auto& D2d = this->Ds_onebase[is][iat1][std::make_pair(iat2, cell)];
for (int iw1 = 0;iw1 < ucell.atoms[it1].nw;++iw1)
for (int iw2 = 0;iw2 < ucell.atoms[it2].nw;++iw2)
D2d(iw1, iw2) = frac * std::conj(this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1))) * this->psi_ks_full(ik, iv, ucell.itiaiw2iwt(it2, ia2, iw2));
D2d(iw1, iw2) = frac * std::conj(this->psi_ks_full(ik, io, ucell.itiaiw2iwt(it1, ia1, iw1))) * this->psi_ks_full(ik, nocc + iv, ucell.itiaiw2iwt(it2, ia2, iw2));
}
}
}
Expand All @@ -83,6 +84,19 @@ namespace LR
psi_out_bfirst.fix_b(ib);
// suppose Cs，Vs, have already been calculated in the ion-step of ground state,
// DM_trans(k) and DM_trans(R) has already been calculated from psi_in in OperatorLRHxc::act
// but int RI_benchmark, DM_trans(k) should be first calculated here
if (cal_dm_trans)
{
#ifdef __MPI
std::vector<container::Tensor> dm_trans_2d = cal_dm_trans_pblas(psi_in_bfirst, *pX, *psi_ks, *pc, naos, nocc, nvirt, *pmat);
if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos, *pmat);
#else
std::vector<container::Tensor> dm_trans_2d = cal_dm_trans_blas(psi_in_bfirst, *psi_ks, nocc, nvirt);
if (this->tdm_sym) for (auto& t : dm_trans_2d) LR_Util::matsym(t.data<T>(), naos);
#endif
// tensor to vector, then set DMK
for (int ik = 0;ik < nk;++ik) { this->DM_trans[ib]->set_DMK_pointer(ik, dm_trans_2d[ik].data<T>()); }
}

// 1. set_Ds (once)
// convert to vector<T*> for the interface of RI_2D_Comm::split_m2D_ktoR (interface will be unified to ct::Tensor)
Expand All @@ -92,17 +106,18 @@ namespace LR
for (int ik = 0;ik < nk;++ik) {DMk_trans_pointer[ik] = &DMk_trans_vector[ik];}
// if multi-k, DM_trans(TR=double) -> Ds_trans(TR=T=complex<double>)
std::vector<std::map<TA, std::map<TAC, RI::Tensor<T>>>> Ds_trans =
RI_2D_Comm::split_m2D_ktoR<T>(this->kv, DMk_trans_pointer, *this->pmat, this->nspin_solve);
RI_2D_Comm::split_m2D_ktoR<T>(this->kv, DMk_trans_pointer, *this->pmat, this->nspin_solve); //0.5 will be multiplied

// 2. cal_Hs
auto lri = this->exx_lri.lock();
for (int ik = 0;ik < nk;++ik)
for (int is = 0;is < nspin_solve;++is)
{
lri->exx_lri.set_Ds(std::move(Ds_trans[ik]), lri->info.dm_threshold);
// LR_Util::print_CV(Ds_trans[is], "Ds_trans in OperatorLREXX", 1e-10);
lri->exx_lri.set_Ds(std::move(Ds_trans[is]), lri->info.dm_threshold);
lri->exx_lri.cal_Hs();
lri->Hexxs[ik] = RI::Communicate_Tensors_Map_Judge::comm_map2_first(
lri->mpi_comm, std::move(lri->exx_lri.Hs), std::get<0>(judge[ik]), std::get<1>(judge[ik]));
lri->post_process_Hexx(lri->Hexxs[ik]);
lri->Hexxs[is] = RI::Communicate_Tensors_Map_Judge::comm_map2_first(
lri->mpi_comm, std::move(lri->exx_lri.Hs), std::get<0>(judge[is]), std::get<1>(judge[is]));
lri->post_process_Hexx(lri->Hexxs[is]);
}

// 3. set [AX]_iak = DM_onbase * Hexxs for each occ-virt pair and each k-point
Expand All @@ -117,12 +132,14 @@ namespace LR
for (int is = 0;is < this->nspin_solve;++is)
{
this->cal_DM_onebase(this->pX->local2global_col(io), this->pX->local2global_row(iv), ik, is); //set Ds_onebase
psi_out_bfirst(ik, io * this->pX->get_row_size() + iv) -= 0.5 * //minus for exchange, 0.5 for spin
// LR_Util::print_CV(Ds_onebase[is], "Ds_onebase of occ " + std::to_string(io) + ", virtual " + std::to_string(iv) + " in OperatorLREXX", 1e-10);
psi_out_bfirst(ik, io * this->pX->get_row_size() + iv) += 2 * //minus for exchange(but here plus is right, why?), 2 for Hartree to Ry
alpha * lri->exx_lri.post_2D.cal_energy(this->Ds_onebase[is], lri->Hexxs[is]);
}
}
}
}

}
}
template class OperatorLREXX<double>;
Expand Down
7 changes: 5 additions & 2 deletions source/module_lr/operator_casida/operator_lr_exx.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,11 @@ namespace LR
Parallel_2D* pX_in,
Parallel_2D* pc_in,
Parallel_Orbitals* pmat_in,
const double& alpha = 1.0)
const double& alpha = 1.0,
const bool& cal_dm_trans = false)
: nspin(nspin), naos(naos), nocc(nocc), nvirt(nvirt),
psi_ks(psi_ks_in), DM_trans(DM_trans_in), exx_lri(exx_lri_in), kv(kv_in),
pX(pX_in), pc(pc_in), pmat(pmat_in), ucell(ucell_in), alpha(alpha)
pX(pX_in), pc(pc_in), pmat(pmat_in), ucell(ucell_in), alpha(alpha), cal_dm_trans(cal_dm_trans)
{
ModuleBase::TITLE("OperatorLREXX", "OperatorLREXX");
this->cal_type = hamilt::calculation_type::lcao_exx;
Expand Down Expand Up @@ -64,6 +65,8 @@ namespace LR
const int& nocc;
const int& nvirt;
const double& alpha;
const bool cal_dm_trans = false;
const bool tdm_sym = false; ///< whether transition density matrix is symmetric
const K_Vectors& kv;
/// ground state wavefunction
const psi::Psi<T>* psi_ks = nullptr;
Expand Down

0 comments on commit 28fa245

Please sign in to comment.