Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor: Unify the interfaces of Pulay terms of force and stress #5130

Merged
merged 8 commits into from
Oct 16, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 3 additions & 6 deletions source/Makefile.Objects
Original file line number Diff line number Diff line change
Expand Up @@ -568,12 +568,9 @@ OBJS_LCAO=evolve_elec.o\
FORCE_STRESS.o\
FORCE_gamma.o\
FORCE_k.o\
fvl_dphi_gamma.o\
fvl_dphi_k.o\
fedm_gamma.o\
fedm_k.o\
ftvnl_dphi_gamma.o\
ftvnl_dphi_k.o\
stress_tools.o\
edm.o\
pulay_force_stress_center2.o\
fvnl_dbeta_gamma.o\
fvnl_dbeta_k.o\
grid_init.o\
Expand Down
9 changes: 3 additions & 6 deletions source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -15,15 +15,12 @@ if(ENABLE_LCAO)
operator_lcao/td_nonlocal_lcao.cpp
operator_lcao/sc_lambda_lcao.cpp
operator_lcao/dftu_lcao.cpp
pulay_force_stress_center2.cpp
FORCE_STRESS.cpp
FORCE_gamma.cpp
FORCE_k.cpp
fvl_dphi_gamma.cpp
fvl_dphi_k.cpp
fedm_gamma.cpp
fedm_k.cpp
ftvnl_dphi_gamma.cpp
ftvnl_dphi_k.cpp
stress_tools.cpp
edm.cpp
fvnl_dbeta_gamma.cpp
fvnl_dbeta_k.cpp
grid_init.cpp
Expand Down
18 changes: 11 additions & 7 deletions 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,12 +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;
};

// this namespace used to store global function for some stress operation
namespace StressTools
{
// set upper matrix to whole matrix
void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix);
} // namespace StressTools
#endif
36 changes: 12 additions & 24 deletions source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "module_elecstate/elecstate_lcao.h"
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
#include "module_io/write_HS.h"
#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"

template <>
void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
Expand Down Expand Up @@ -215,10 +216,17 @@ void Force_LCAO<double>::ftable(const bool isforce,
// allocate DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z
this->allocate(pv, fsr, two_center_bundle, orb);

const double* dSx[3] = { fsr.DSloc_x, fsr.DSloc_y, fsr.DSloc_z };
const double* dSxy[6] = { fsr.DSloc_11, fsr.DSloc_12, fsr.DSloc_13, fsr.DSloc_22, fsr.DSloc_23, fsr.DSloc_33 };
// calculate the force related to 'energy density matrix'.
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap);
PulayForceStress::cal_pulay_fs(foverlap, soverlap,
this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
ucell, pv, dSx, dSxy, isforce, isstress);

this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);
const double* dHx[3] = { fsr.DHloc_fixed_x, fsr.DHloc_fixed_y, fsr.DHloc_fixed_z };
const double* dHxy[6] = { fsr.DHloc_fixed_11, fsr.DHloc_fixed_12, fsr.DHloc_fixed_13, fsr.DHloc_fixed_22, fsr.DHloc_fixed_23, fsr.DHloc_fixed_33 };
//tvnl_dphi
PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);

this->cal_fvnl_dbeta(dm,
pv,
Expand All @@ -231,9 +239,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
fvnl_dbeta,
svnl_dbeta);

this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi);
// vl_dphi
PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);

// caoyu add for DeePKS
#ifdef __DEEPKS
if (PARAM.inp.deepks_scf)
{
Expand Down Expand Up @@ -314,23 +322,3 @@ void Force_LCAO<double>::ftable(const bool isforce,
ModuleBase::timer::tick("Force_LCAO", "ftable");
return;
}

namespace StressTools
{
void stress_fill(const double& lat0_, const double& omega_, ModuleBase::matrix& stress_matrix)
{
assert(omega_ > 0.0);
double weight = lat0_ / omega_;
for (int i = 0; i < 3; ++i)
{
for (int j = 0; j < 3; ++j)
{
if (j > i)
{
stress_matrix(j, i) = stress_matrix(i, j);
}
stress_matrix(i, j) *= weight;
}
}
}
} // namespace StressTools
14 changes: 11 additions & 3 deletions source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h"
#include "module_hamilt_pw/hamilt_pwdft/global.h"
#include "module_io/write_HS.h"
#include "module_hamilt_lcao/hamilt_lcaodft/pulay_force_stress.h"

#include <map>
#include <unordered_map>
Expand Down Expand Up @@ -312,14 +313,21 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
kv->get_nks(),
kv->kvec_d);

const double* dSx[3] = { fsr.DSloc_Rx, fsr.DSloc_Ry, fsr.DSloc_Rz };
// 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);
PulayForceStress::cal_pulay_fs(foverlap, soverlap,
this->cal_edm(pelec, *psi, *dm, *kv, pv, PARAM.inp.nspin, PARAM.inp.nbands, ucell, *ra),
ucell, pv, dSx, fsr.DH_r, isforce, isstress, ra, -1.0, 1.0);

this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
const double* dHx[3] = { fsr.DHloc_fixedR_x, fsr.DHloc_fixedR_y, fsr.DHloc_fixedR_z }; // T+Vnl
const double* dHxy[6] = { fsr.stvnl11, fsr.stvnl12, fsr.stvnl13, fsr.stvnl22, fsr.stvnl23, fsr.stvnl33 }; //T
// tvnl_dphi
PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress, ra, 1.0, -1.0);

// doing on the real space grid.
this->cal_fvl_dphi(isforce, isstress, pelec->pot, gint, fvl_dphi, svl_dphi);
// vl_dphi
PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);

this->cal_fvnl_dbeta(dm,
pv,
Expand Down
103 changes: 103 additions & 0 deletions source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
#include "FORCE.h"
#include "module_elecstate/module_dm/cal_dm_psi.h"
#include "module_base/memory.h"
#include "module_parameter/parameter.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 (PARAM.inp.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;
}
112 changes: 0 additions & 112 deletions source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp

This file was deleted.

Loading
Loading