Skip to content
2 changes: 1 addition & 1 deletion source/source_io/output_mat_sparse.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ void output_mat_sparse(const bool& out_mat_hsR,
//! generate a file containing the Hamiltonian and S(overlap) matrices
if (out_mat_hsR)
{
output_HSR(ucell, istep, v_eff, pv, HS_Arrays, grid, kv, *p_dftu, p_ham);
output_HSR(ucell, istep, pv, HS_Arrays, grid, kv, *p_dftu, p_ham);
}

//! generate a file containing the kinetic energy matrix
Expand Down
3 changes: 0 additions & 3 deletions source/source_io/write_HS_R.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@
template <typename TK>
void ModuleIO::output_HSR(const UnitCell& ucell,
const int& istep,
const ModuleBase::matrix& v_eff,
const Parallel_Orbitals& pv,
LCAO_HS_Arrays& HS_Arrays,
const Grid_Driver& grid, // mohan add 2024-04-06
Expand Down Expand Up @@ -304,7 +303,6 @@ void ModuleIO::output_TR(const int istep,
template void ModuleIO::output_HSR<double>(
const UnitCell& ucell,
const int& istep,
const ModuleBase::matrix& v_eff,
const Parallel_Orbitals& pv,
LCAO_HS_Arrays& HS_Arrays,
const Grid_Driver& grid,
Expand All @@ -324,7 +322,6 @@ template void ModuleIO::output_HSR<double>(
template void ModuleIO::output_HSR<std::complex<double>>(
const UnitCell& ucell,
const int& istep,
const ModuleBase::matrix& v_eff,
const Parallel_Orbitals& pv,
LCAO_HS_Arrays& HS_Arrays,
const Grid_Driver& grid,
Expand Down
1 change: 0 additions & 1 deletion source/source_io/write_HS_R.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ using TAC = std::pair<int, std::array<int, 3>>;
template <typename TK>
void output_HSR(const UnitCell& ucell,
const int& istep,
const ModuleBase::matrix& v_eff,
const Parallel_Orbitals& pv,
LCAO_HS_Arrays& HS_Arrays,
const Grid_Driver& grid, // mohan add 2024-04-06
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/FORCE_gamma.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -231,10 +231,10 @@ void Force_LCAO<double>::ftable(const bool isforce,
gd,
*this->ParaV,
nks,
deepks.ld.deepks_param,
kv->kvec_d,
deepks.ld.phialpha,
deepks.ld.gedm,
deepks.ld.inl_index,
fvnl_dalpha,
isstress,
svnl_dalpha);
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/FORCE_k.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -257,10 +257,10 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
gd,
pv,
kv->get_nks(),
deepks.ld.deepks_param,
kv->kvec_d,
deepks.ld.phialpha,
deepks.ld.gedm,
deepks.ld.inl_index,
fvnl_dalpha,
isstress,
svnl_dalpha);
Expand Down
88 changes: 39 additions & 49 deletions source/source_lcao/LCAO_hamilt.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@
#include "source_base/abfs-vector3_order.h"
#include "source_base/global_variable.h"
#include "source_base/timer.h"
#include "source_lcao/spar_exx.h"
#include "source_lcao/module_ri/RI_2D_Comm.h"
#include "source_lcao/spar_exx.h"

#include <RI/global/Global_Func-2.h>
#include <RI/ri/Cell_Nearest.h>
Expand All @@ -21,28 +21,27 @@
// Peize Lin add 2022.09.13

template <typename Tdata>
void sparse_format::cal_HR_exx(const UnitCell& ucell,
void sparse_format::cal_HR_exx(
const UnitCell& ucell,
const Parallel_Orbitals& pv,
LCAO_HS_Arrays& HS_Arrays,
const int& current_spin,
const double& sparse_threshold,
const int (&nmp)[3],
const std::vector<std::map<int,
std::map<std::pair<int, std::array<int, 3>>,
RI::Tensor<Tdata>>>>& Hexxs) {
const std::vector<std::map<int, std::map<std::pair<int, std::array<int, 3>>, RI::Tensor<Tdata>>>>& Hexxs)
{
ModuleBase::TITLE("sparse_format", "cal_HR_exx");
ModuleBase::timer::tick("sparse_format", "cal_HR_exx");

const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha;

std::map<int, std::array<double, 3>> atoms_pos;
for (int iat = 0; iat < ucell.nat; ++iat) {
atoms_pos[iat] = RI_Util::Vector3_to_array3(
ucell.atoms[ucell.iat2it[iat]]
.tau[ucell.iat2ia[iat]]);
for (int iat = 0; iat < ucell.nat; ++iat)
{
atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]);
}
const std::array<std::array<double, 3>, 3> latvec
= {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here,
= {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here,
RI_Util::Vector3_to_array3(ucell.a2),
RI_Util::Vector3_to_array3(ucell.a3)};

Expand All @@ -51,92 +50,83 @@ void sparse_format::cal_HR_exx(const UnitCell& ucell,
RI::Cell_Nearest<int, int, 3, double, 3> cell_nearest;
cell_nearest.init(atoms_pos, latvec, Rs_period);

const std::vector<int> is_list = (PARAM.inp.nspin != 4)
? std::vector<int>{current_spin}
: std::vector<int>{0, 1, 2, 3};
const std::vector<int> is_list
= (PARAM.inp.nspin != 4) ? std::vector<int>{current_spin} : std::vector<int>{0, 1, 2, 3};

for (const int is: is_list)
for (const int is: is_list)
{
int is0_b = 0;
int is1_b = 0;
std::tie(is0_b, is1_b) = RI_2D_Comm::split_is_block(is);

if (Hexxs.empty())
if (Hexxs.empty())
{
break;
}

for (const auto& HexxA: Hexxs[is])
for (const auto& HexxA: Hexxs[is])
{
const int iat0 = HexxA.first;
for (const auto& HexxB: HexxA.second)
for (const auto& HexxB: HexxA.second)
{
const int iat1 = HexxB.first.first;

const Abfs::Vector3_Order<int> R = RI_Util::array3_to_Vector3(
cell_nearest.get_cell_nearest_discrete(iat0,
iat1,
HexxB.first.second));
cell_nearest.get_cell_nearest_discrete(iat0, iat1, HexxB.first.second));

HS_Arrays.all_R_coor.insert(R);

const RI::Tensor<Tdata>& Hexx = HexxB.second;

for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0)
for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0)
{
const int iwt0 = RI_2D_Comm::get_iwt(ucell,iat0, iw0, is0_b);
const int iwt0 = RI_2D_Comm::get_iwt(ucell, iat0, iw0, is0_b);
const int iwt0_local = pv.global2local_row(iwt0);

if (iwt0_local < 0)
if (iwt0_local < 0)
{
continue;
}

for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1)
for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1)
{
const int iwt1 = RI_2D_Comm::get_iwt(ucell,iat1, iw1, is1_b);
const int iwt1 = RI_2D_Comm::get_iwt(ucell, iat1, iw1, is1_b);
const int iwt1_local = pv.global2local_col(iwt1);

if (iwt1_local < 0)
if (iwt1_local < 0)
{
continue;
}

if (std::abs(Hexx(iw0, iw1)) > sparse_threshold)
if (std::abs(Hexx(iw0, iw1)) > sparse_threshold)
{
if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2)
if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2)
{
auto& HR_sparse_ptr
= HS_Arrays
.HR_sparse[current_spin][R][iwt0];
auto& HR_sparse_ptr = HS_Arrays.HR_sparse[current_spin][R][iwt0];
double& HR_sparse = HR_sparse_ptr[iwt1];
HR_sparse += RI::Global_Func::convert<double>(
frac * Hexx(iw0, iw1));
if (std::abs(HR_sparse) <= sparse_threshold)
HR_sparse += RI::Global_Func::convert<double>(frac * Hexx(iw0, iw1));
if (std::abs(HR_sparse) <= sparse_threshold)
{
HR_sparse_ptr.erase(iwt1);
}
}
else if (PARAM.inp.nspin == 4)
}
else if (PARAM.inp.nspin == 4)
{
auto& HR_sparse_ptr
= HS_Arrays.HR_soc_sparse[R][iwt0];

std::complex<double>& HR_sparse
= HR_sparse_ptr[iwt1];

HR_sparse += RI::Global_Func::convert<
std::complex<double>>(frac * Hexx(iw0, iw1));

if (std::abs(HR_sparse) <= sparse_threshold)
auto& HR_sparse_ptr = HS_Arrays.HR_soc_sparse[R][iwt0];

std::complex<double>& HR_sparse = HR_sparse_ptr[iwt1];

HR_sparse += RI::Global_Func::convert<std::complex<double>>(frac * Hexx(iw0, iw1));

if (std::abs(HR_sparse) <= sparse_threshold)
{
HR_sparse_ptr.erase(iwt1);
}
}
else
}
else
{
throw std::invalid_argument(std::string(__FILE__) + " line "
+ std::to_string(__LINE__));
+ std::to_string(__LINE__));
}
}
}
Expand Down
Loading
Loading