Skip to content

Commit

Permalink
Merge branch 'develop' into rho-pcie
Browse files Browse the repository at this point in the history
  • Loading branch information
@dzzz2001
dzzz2001 authored Jul 3, 2024
2 parents 2135a03 + 3726d73 commit 57f5e4c
Show file tree
Hide file tree
Showing 21 changed files with 1,833 additions and 1,795 deletions.
4 changes: 1 addition & 3 deletions source/Makefile.Objects
View file
Original file line number Diff line number Diff line change
Expand Up @@ -477,16 +477,14 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
nscf_fermi_surf.o\
istate_charge.o\
istate_envelope.o\
read_dm.o\
io_dmk.o\
unk_overlap_lcao.o\
read_wfc_nao.o\
read_wfc_lcao.o\
write_wfc_nao.o\
write_HS_sparse.o\
single_R_io.o\
write_HS_R.o\
write_dm.o\
output_dm.o\
write_dmr.o\
sparse_matrix.o\
output_mulliken.o\
Expand Down
172 changes: 77 additions & 95 deletions source/module_elecstate/elecstate_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,72 +10,78 @@

#include <vector>

namespace elecstate
{
namespace elecstate {

// multi-k case
template <>
void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<double>>& psi)
{
void ElecStateLCAO<std::complex<double>>::psiToRho(
const psi::Psi<std::complex<double>>& psi) {
ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");

this->calculate_weights();

// the calculations of dm, and dm -> rho are, technically, two separate functionalities, as we cannot
// rule out the possibility that we may have a dm from other sources, such as read from file.
// However, since we are not separating them now, I opt to add a flag to control how dm is obtained as of now
if (!GlobalV::dm_to_rho)
{
// the calculations of dm, and dm -> rho are, technically, two separate
// functionalities, as we cannot rule out the possibility that we may have a
// dm from other sources, such as read from file. However, since we are not
// separating them now, I opt to add a flag to control how dm is obtained as
// of now
if (!GlobalV::dm_to_rho) {
this->calEBand();

ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");

// this part for calculating DMK in 2d-block format, not used for charge now
// this part for calculating DMK in 2d-block format, not used for charge
// now
// psi::Psi<std::complex<double>> dm_k_2d();

if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack"
|| GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp"
if (GlobalV::KS_SOLVER == "genelpa"
|| GlobalV::KS_SOLVER == "scalapack_gvx"
|| GlobalV::KS_SOLVER == "lapack"
|| GlobalV::KS_SOLVER == "cusolver"
|| GlobalV::KS_SOLVER == "cusolvermp"
|| GlobalV::KS_SOLVER == "cg_in_lcao") // Peize Lin test 2019-05-15
{
// cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_k);
elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM));
elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
this->wg,
psi,
*(this->DM));
this->DM->cal_DMR();

// interface for RI-related calculation, which needs loc.dm_k
#ifdef __EXX
if (GlobalC::exx_info.info_global.cal_exx)
{
if (GlobalC::exx_info.info_global.cal_exx) {
const K_Vectors* kv = this->DM->get_kv_pointer();
this->loc->dm_k.resize(kv->get_nks());
for (int ik = 0; ik < kv->get_nks(); ++ik)
{
for (int ik = 0; ik < kv->get_nks(); ++ik) {
this->loc->set_dm_k(ik, this->DM->get_DMK_pointer(ik));
}
}
#endif
}
}

for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],
this->charge->nrxx); // mohan 2009-11-10
}

//------------------------------------------------------------
// calculate the charge density on real space grid.
//------------------------------------------------------------

ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
this->gint_k->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
this->gint_k->transfer_DM2DtoGrid(
this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho);
this->gint_k->cal_gint(&inout);

if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
{
for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
if (XC_Functional::get_func_type() == 3
|| XC_Functional::get_func_type() == 5) {
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is],
this->charge->nrxx);
}
Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
this->gint_k->cal_gint(&inout1);
Expand All @@ -89,64 +95,50 @@ void ElecStateLCAO<std::complex<double>>::psiToRho(const psi::Psi<std::complex<d

// Gamma_only case
template <>
void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
{
void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi) {
ModuleBase::TITLE("ElecStateLCAO", "psiToRho");
ModuleBase::timer::tick("ElecStateLCAO", "psiToRho");

this->calculate_weights();
this->calEBand();

if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx" || GlobalV::KS_SOLVER == "lapack"
|| GlobalV::KS_SOLVER == "cusolver" || GlobalV::KS_SOLVER == "cusolvermp" || GlobalV::KS_SOLVER == "cg_in_lcao")
{
if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx"
|| GlobalV::KS_SOLVER == "lapack" || GlobalV::KS_SOLVER == "cusolver"
|| GlobalV::KS_SOLVER == "cusolvermp"
|| GlobalV::KS_SOLVER == "cg_in_lcao") {
ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");

// get DMK in 2d-block format
// cal_dm(this->loc->ParaV, this->wg, psi, this->loc->dm_gamma);
elecstate::cal_dm_psi(this->DM->get_paraV_pointer(), this->wg, psi, *(this->DM));
elecstate::cal_dm_psi(this->DM->get_paraV_pointer(),
this->wg,
psi,
*(this->DM));
this->DM->cal_DMR();
if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready
{
this->loc->dm_gamma.resize(GlobalV::NSPIN);
for (int is = 0; is < GlobalV::NSPIN; ++is)
{
this->loc->set_dm_gamma(is, this->DM->get_DMK_pointer(is));
}
}
ModuleBase::timer::tick("ElecStateLCAO", "cal_dm_2d");
for (int ik = 0; ik < psi.get_nk(); ++ik)
{
// for gamma_only case, no convertion occured, just for print.
// old 2D-to-Grid conversion has been replaced by new Gint Refactor 2023/09/25
if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready
{
this->loc->cal_dk_gamma_from_2D_pub();
}
}
}

for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],
this->charge->nrxx); // mohan 2009-11-10
}

//------------------------------------------------------------
// calculate the charge density on real space grid.
//------------------------------------------------------------
ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");

this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
this->gint_gamma->transfer_DM2DtoGrid(
this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint

Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho);

this->gint_gamma->cal_gint(&inout);

if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
{
for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is], this->charge->nrxx);
if (XC_Functional::get_func_type() == 3
|| XC_Functional::get_func_type() == 5) {
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[is],
this->charge->nrxx);
}
Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
this->gint_gamma->cal_gint(&inout1);
Expand All @@ -159,70 +151,59 @@ void ElecStateLCAO<double>::psiToRho(const psi::Psi<double>& psi)
}

template <typename TK>
void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv, const Parallel_Orbitals* paraV, const int nspin)
{
void ElecStateLCAO<TK>::init_DM(const K_Vectors* kv,
const Parallel_Orbitals* paraV,
const int nspin) {
this->DM = new DensityMatrix<TK, double>(kv, paraV, nspin);
}

template <>
double ElecStateLCAO<double>::get_spin_constrain_energy()
{
SpinConstrain<double, base_device::DEVICE_CPU>& sc = SpinConstrain<double>::getScInstance();
double ElecStateLCAO<double>::get_spin_constrain_energy() {
SpinConstrain<double, base_device::DEVICE_CPU>& sc
= SpinConstrain<double>::getScInstance();
return sc.cal_escon();
}

template <>
double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy()
{
double ElecStateLCAO<std::complex<double>>::get_spin_constrain_energy() {
SpinConstrain<std::complex<double>, base_device::DEVICE_CPU>& sc
= SpinConstrain<std::complex<double>>::getScInstance();
return sc.cal_escon();
}

#ifdef __PEXSI
template <>
void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM, std::vector<double*> pexsi_EDM)
{
void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM,
std::vector<double*> pexsi_EDM) {
ModuleBase::timer::tick("ElecStateLCAO", "dmToRho");

int nspin = GlobalV::NSPIN;
if (GlobalV::NSPIN == 4)
{
if (GlobalV::NSPIN == 4) {
nspin = 1;
}

// old 2D-to-Grid conversion has been replaced by new Gint Refactor 2023/09/25
if (this->loc->out_dm) // keep interface for old Output_DM until new one is ready
{
for (int is = 0; is < nspin; ++is)
{
this->loc->set_dm_gamma(is, pexsi_DM[is]);
}
this->loc->cal_dk_gamma_from_2D_pub();
}

this->get_DM()->pexsi_EDM = pexsi_EDM;

for (int is = 0; is < nspin; is++)
{
for (int is = 0; is < nspin; is++) {
this->DM->set_DMK_pointer(is, pexsi_DM[is]);
}
DM->cal_DMR();

for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is], this->charge->nrxx); // mohan 2009-11-10
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->rho[is],
this->charge->nrxx); // mohan 2009-11-10
}

ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
this->gint_gamma->transfer_DM2DtoGrid(this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
this->gint_gamma->transfer_DM2DtoGrid(
this->DM->get_DMR_vector()); // transfer DM2D to DM_grid in gint
Gint_inout inout(this->charge->rho, Gint_Tools::job_type::rho);
this->gint_gamma->cal_gint(&inout);
if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
{
for (int is = 0; is < GlobalV::NSPIN; is++)
{
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0], this->charge->nrxx);
if (XC_Functional::get_func_type() == 3
|| XC_Functional::get_func_type() == 5) {
for (int is = 0; is < GlobalV::NSPIN; is++) {
ModuleBase::GlobalFunc::ZEROS(this->charge->kin_r[0],
this->charge->nrxx);
}
Gint_inout inout1(this->charge->kin_r, Gint_Tools::job_type::tau);
this->gint_gamma->cal_gint(&inout1);
Expand All @@ -235,10 +216,11 @@ void ElecStateLCAO<double>::dmToRho(std::vector<double*> pexsi_DM, std::vector<d
}

template <>
void ElecStateLCAO<std::complex<double>>::dmToRho(std::vector<std::complex<double>*> pexsi_DM,
std::vector<std::complex<double>*> pexsi_EDM)
{
ModuleBase::WARNING_QUIT("ElecStateLCAO", "pexsi is not completed for multi-k case");
void ElecStateLCAO<std::complex<double>>::dmToRho(
std::vector<std::complex<double>*> pexsi_DM,
std::vector<std::complex<double>*> pexsi_EDM) {
ModuleBase::WARNING_QUIT("ElecStateLCAO",
"pexsi is not completed for multi-k case");
}

#endif
Expand Down
Loading

0 comments on commit 57f5e4c

Please sign in to comment.