Skip to content

Commit

Permalink
extract edm
Browse files Browse the repository at this point in the history
  • Loading branch information
maki49 committed Sep 17, 2024
1 parent efb0990 commit 86e0291
Show file tree
Hide file tree
Showing 8 changed files with 125 additions and 84 deletions.
1 change: 1 addition & 0 deletions source/Makefile.Objects
Original file line number Diff line number Diff line change
Expand Up @@ -560,6 +560,7 @@ OBJS_LCAO=evolve_elec.o\
FORCE_gamma.o\
FORCE_k.o\
stress_tools.o\
edm.o\
fedm_gamma.o\
fedm_k.o\
ftvnl_dphi_gamma.o\
Expand Down
1 change: 1 addition & 0 deletions source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ if(ENABLE_LCAO)
FORCE_gamma.cpp
FORCE_k.cpp
stress_tools.cpp
edm.cpp
fedm_gamma.cpp
fedm_k.cpp
ftvnl_dphi_gamma.cpp
Expand Down
12 changes: 11 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ class Force_LCAO
const bool isstress,
ForceStressArrays& fsr,
const UnitCell& ucell,
const elecstate::DensityMatrix<T, double>* dm,
const elecstate::DensityMatrix<T, double>& dm,
const psi::Psi<T>* psi,
const Parallel_Orbitals& pv,
const elecstate::ElecState* pelec,
Expand Down Expand Up @@ -136,6 +136,16 @@ class Force_LCAO
typename TGint<T>::type& gint,
ModuleBase::matrix& fvl_dphi,
ModuleBase::matrix& svl_dphi);

elecstate::DensityMatrix<T, double> cal_edm(const elecstate::ElecState* pelec,
const psi::Psi<T>& psi,
const elecstate::DensityMatrix<T, double>& dm,
const K_Vectors& kv,
const Parallel_Orbitals& pv,
const int& nspin,
const int& nbands,
const UnitCell& ucell,
Record_adj& ra) const;
};

#endif
4 changes: 3 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
this->allocate(pv, fsr, two_center_bundle, orb);

// calculate the force related to 'energy density matrix'.
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap);
this->cal_fedm(isforce, isstress, fsr, ucell,
this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra),
psi, pv, pelec, foverlap, soverlap);

this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);

Expand Down
4 changes: 3 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -312,7 +312,9 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,

// calculate the energy density matrix
// and the force related to overlap matrix and energy density matrix.
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap, kv, ra);
this->cal_fedm(isforce, isstress, fsr, ucell,
this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra),
psi, pv, pelec, foverlap, soverlap, kv, ra);

this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);

Expand Down
102 changes: 102 additions & 0 deletions source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
#include "FORCE.h"
#include "module_elecstate/module_dm/cal_dm_psi.h"
#include "module_base/memory.h"
template<>
elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecstate::ElecState* pelec,
const psi::Psi<double>& psi,
const elecstate::DensityMatrix<double, double>& dm,
const K_Vectors& kv,
const Parallel_Orbitals& pv,
const int& nspin,
const int& nbands,
const UnitCell& ucell,
Record_adj& ra) const
{
ModuleBase::matrix wg_ekb;
wg_ekb.create(nspin, nbands);

for(int is=0; is<nspin; is++)
{
for(int ib=0; ib<nbands; ib++)
{
wg_ekb(is,ib) = pelec->wg(is,ib) * pelec->ekb(is, ib);
}
}

// construct a DensityMatrix for Gamma-Only
elecstate::DensityMatrix<double, double> edm(&pv, nspin);

#ifdef __PEXSI
if (GlobalV::KS_SOLVER == "pexsi")
{
auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
for (int ik = 0; ik < nspin; ik++)
{
edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
}

}
else
#endif
{
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
}
return edm;
}

template<>
elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(const elecstate::ElecState* pelec,
const psi::Psi<std::complex<double>>& psi,
const elecstate::DensityMatrix<std::complex<double>, double>& dm,
const K_Vectors& kv,
const Parallel_Orbitals& pv,
const int& nspin,
const int& nbands,
const UnitCell& ucell,
Record_adj& ra) const
{

// construct a DensityMatrix object
elecstate::DensityMatrix<std::complex<double>, double> edm(&kv, &pv, nspin);

//--------------------------------------------
// calculate the energy density matrix here.
//--------------------------------------------

ModuleBase::matrix wg_ekb;
wg_ekb.create(kv.get_nks(), nbands);
ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv.get_nks() * nbands);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 1024)
#endif
for (int ik = 0; ik < kv.get_nks(); ik++)
{
for (int ib = 0; ib < nbands; ib++)
{
wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib);
}
}

// use the original formula (Hamiltonian matrix) to calculate energy density matrix
if (dm.EDMK.size())
{
#ifdef _OPENMP
#pragma omp parallel for schedule(static)
#endif
for (int ik = 0; ik < kv.get_nks(); ++ik)
{
edm.set_DMK_pointer(ik, dm.EDMK[ik].c);
}
}
else
{
// cal_dm_psi
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
}

// cal_dm_2d
edm.init_DMR(ra, &ucell);
edm.cal_DMR();
edm.sum_DMR_spin();
return edm;
}
35 changes: 2 additions & 33 deletions source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ void Force_LCAO<double>::cal_fedm(
const bool isstress,
ForceStressArrays &fsr, // mohan add 2024-06-16
const UnitCell& ucell,
const elecstate::DensityMatrix<double, double>* dm,
const elecstate::DensityMatrix<double, double>& dm,
const psi::Psi<double>* psi,
const Parallel_Orbitals& pv,
const elecstate::ElecState *pelec,
Expand All @@ -32,37 +32,6 @@ void Force_LCAO<double>::cal_fedm(
const int nbands = GlobalV::NBANDS;
const int nlocal = GlobalV::NLOCAL;

ModuleBase::matrix wg_ekb;
wg_ekb.create(nspin, nbands);

for(int is=0; is<nspin; is++)
{
for(int ib=0; ib<nbands; ib++)
{
wg_ekb(is,ib) = pelec->wg(is,ib) * pelec->ekb(is, ib);
}
}

// construct a DensityMatrix for Gamma-Only
elecstate::DensityMatrix<double, double> edm(&pv, nspin);

#ifdef __PEXSI
if (GlobalV::KS_SOLVER == "pexsi")
{
auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
for (int ik = 0; ik < nspin; ik++)
{
edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
}

}
else
#endif
{
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm);
}


for(int i=0; i<nlocal; i++)
{
const int iat = ucell.iwt2iat[i];
Expand All @@ -78,7 +47,7 @@ void Force_LCAO<double>::cal_fedm(
for(int is=0; is<nspin; ++is)
{
//sum += edm_gamma[is](nu, mu);
sum += edm.get_DMK(is+1, 0, nu, mu);
sum += dm.get_DMK(is+1, 0, nu, mu);
}
sum *= 2.0;

Expand Down
50 changes: 2 additions & 48 deletions source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
#include "module_cell/module_neighbor/sltk_grid_driver.h"
#include "module_elecstate/cal_dm.h"
#include "module_elecstate/elecstate_lcao.h"
#include "module_elecstate/module_dm/cal_dm_psi.h"
#include "module_hamilt_pw/hamilt_pwdft/global.h"
#include "module_io/write_HS.h"
#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
Expand All @@ -29,7 +28,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
const bool isstress,
ForceStressArrays& fsr,
const UnitCell& ucell,
const elecstate::DensityMatrix<std::complex<double>, double>* dm,
const elecstate::DensityMatrix<std::complex<double>, double>& dm,
const psi::Psi<std::complex<double>>* psi,
const Parallel_Orbitals& pv,
const elecstate::ElecState* pelec,
Expand All @@ -42,51 +41,6 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
ModuleBase::timer::tick("Force_LCAO", "cal_fedm");

const int nspin = GlobalV::NSPIN;
const int nbands = GlobalV::NBANDS;

// construct a DensityMatrix object
elecstate::DensityMatrix<std::complex<double>, double> edm(kv, &pv, nspin);

//--------------------------------------------
// calculate the energy density matrix here.
//--------------------------------------------

ModuleBase::matrix wg_ekb;
wg_ekb.create(kv->get_nks(), nbands);
ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv->get_nks() * nbands);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 1024)
#endif
for (int ik = 0; ik < kv->get_nks(); ik++)
{
for (int ib = 0; ib < nbands; ib++)
{
wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib);
}
}
std::vector<ModuleBase::ComplexMatrix> edm_k(kv->get_nks());

// use the original formula (Hamiltonian matrix) to calculate energy density matrix
if (dm->EDMK.size())
{
#ifdef _OPENMP
#pragma omp parallel for schedule(static)
#endif
for (int ik = 0; ik < kv->get_nks(); ++ik)
{
edm.set_DMK_pointer(ik, dm->EDMK[ik].c);
}
}
else
{
// cal_dm_psi
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm);
}

// cal_dm_2d
edm.init_DMR(*ra, &ucell);
edm.cal_DMR();
edm.sum_DMR_spin();
//--------------------------------------------
// summation \sum_{i,j} E(i,j)*dS(i,j)
// BEGIN CALCULATION OF FORCE OF EACH ATOM
Expand Down Expand Up @@ -146,7 +100,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
double Ry = ra->info[iat][cb][1];
double Rz = ra->info[iat][cb][2];
// get BaseMatrix
hamilt::BaseMatrix<double>* tmp_matrix = edm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
hamilt::BaseMatrix<double>* tmp_matrix = dm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
if (tmp_matrix == nullptr)
{
continue;
Expand Down

0 comments on commit 86e0291

Please sign in to comment.