Skip to content

Fix(dftu-pw): fix DFT+U for PW basis (nspin=1/2/4) and occupation mat… #7733

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.

Sign up for GitHub

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

Closed
dyzheng wants to merge 2 commits into from
Closed

Fix(dftu-pw): fix DFT+U for PW basis (nspin=1/2/4) and occupation mat… #7733

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6630a2b
Fix(dftu-pw): fix DFT+U for PW basis (nspin=1/2/4) and occupation mat…
Jul 31, 2026
ff43ca7
Merge branch 'develop' into feat/issue-7726-dftu-pw-nspin-fixes
mohanchen Aug 1, 2026
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
12 changes: 4 additions & 8 deletions source/source_lcao/module_dftu/dftu.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,10 @@ void Plus_U::init(UnitCell& cell,
// it:index of type of atom
for (int it = 0; it < cell.ntype; ++it)
{
if(!has_correlated_orbital(it))
{
continue;
}
for (int ia = 0; ia < cell.atoms[it].na; ia++)
{
// ia:index of atoms of this type
Expand All @@ -140,14 +144,6 @@ void Plus_U::init(UnitCell& cell,
locale[iat].resize(cell.atoms[it].nwl + 1);
locale_save[iat].resize(cell.atoms[it].nwl + 1);

// initialize the arrry iatlnm2iwt[iat][l][n][m]
this->iatlnmipol2iwt[iat].resize(cell.atoms[it].nwl + 1);

if(!has_correlated_orbital(it))
{
continue;
}

const int tlp1_npol = (get_orbital_corr(it)*2+1)*npol;
const int tlp1 = 2 * get_orbital_corr(it) + 1;
const int elem_size = tlp1 * tlp1;
Expand Down
2 changes: 1 addition & 1 deletion source/source_lcao/module_dftu/dftu.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
#ifdef __LCAO
#include "source_basis/module_ao/ORB_read.h"
#include "source_hamilt/hamilt.h"
#include "source_hamilt/module_hcontainer/hcontainer.h"
#include "source_lcao/module_hcontainer/hcontainer.h"
#include "source_estate/module_dm/density_matrix.h"
#include "source_lcao/force_stress_arrays.h" // mohan add 2024-06-15
#endif
Expand Down
4 changes: 2 additions & 2 deletions source/source_lcao/module_dftu/dftu_folding.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@
#include "source_io/module_parameter/parameter.h"
#include "source_cell/module_neighbor/sltk_grid_driver.h"
#include "source_lcao/hamilt_lcao.h"
#include "source_hamilt/module_hcontainer/hcontainer.h"
#include "source_hamilt/module_hcontainer/hcontainer_funcs.h"
#include "source_lcao/module_hcontainer/hcontainer.h"
#include "source_lcao/module_hcontainer/hcontainer_funcs.h"

void Plus_U::fold_dSR_gamma(const UnitCell& ucell,
const Parallel_Orbitals& pv,
Expand Down
51 changes: 29 additions & 22 deletions source/source_lcao/module_dftu/dftu_io.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,10 @@ void Plus_U::output(const UnitCell& ucell,

if (L >= get_orbital_corr(T) && has_correlated_orbital(T))
{
if (L != get_orbital_corr(T))
{
continue;
}
if (L != get_orbital_corr(T))
{
continue;
}

if (!Yukawa)
{
Expand Down Expand Up @@ -97,11 +97,11 @@ void Plus_U::write_occup_m(const UnitCell& ucell,

for (int T = 0; T < ucell.ntype; T++)
{
if (!has_correlated_orbital(T))
{
continue;
}
const int NL = ucell.atoms[T].nwl + 1;
if (!has_correlated_orbital(T))
{
continue;
}
const int NL = ucell.atoms[T].nwl + 1;
const int LC = get_orbital_corr(T);

for (int I = 0; I < ucell.atoms[T].na; I++)
Expand All @@ -110,10 +110,10 @@ void Plus_U::write_occup_m(const UnitCell& ucell,

for (int l = 0; l < NL; l++)
{
if (l != get_orbital_corr(T))
{
continue;
}
if (l != get_orbital_corr(T))
{
continue;
}

const int N = ucell.atoms[T].l_nchi[l];

Expand Down Expand Up @@ -316,7 +316,14 @@ void Plus_U::read_occup_m(const UnitCell& ucell,

for (int l = 0; l < NL; l++)
{
if (l != get_orbital_corr(T))
if (l != get_orbital_corr(T))
{
continue;
}

ifdftu >> word;

if (strcmp("L", word) == 0)
{
continue;
}
Expand Down Expand Up @@ -397,10 +404,10 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell,

for (int T = 0; T < ucell.ntype; T++)
{
if (!has_correlated_orbital(T))
{
continue;
}
if (!has_correlated_orbital(T))
{
continue;
}

for (int I = 0; I < ucell.atoms[T].na; I++)
{
Expand All @@ -409,10 +416,10 @@ void Plus_U::local_occup_bcast(const UnitCell& ucell,

for (int l = 0; l <= ucell.atoms[T].nwl; l++)
{
if (l != get_orbital_corr(T))
{
continue;
}
if (l != get_orbital_corr(T))
{
continue;
}

for (int n = 0; n < ucell.atoms[T].l_nchi[l]; n++)
{
Expand Down
10 changes: 7 additions & 3 deletions source/source_lcao/module_dftu/dftu_occup.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ void Plus_U::copy_locale(const UnitCell& ucell)
if (Plus_U::nspin == 4)
{
locale_save[iat][target_l][0][0] = locale[iat][target_l][0][0];
// nspin=4 locale matrix already contains all spin components interleaved
if(this->uom_save.size() != 0)
{
const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc;
Expand All @@ -43,6 +44,7 @@ void Plus_U::copy_locale(const UnitCell& ucell)
{
locale_save[iat][target_l][0][0] = locale[iat][target_l][0][0];
locale_save[iat][target_l][0][1] = locale[iat][target_l][0][1];
// save locale matrix for spin=0,1 to uom_save
if(this->uom_save.size() != 0)
{
const int size = locale[iat][target_l][0][0].nr * locale[iat][target_l][0][0].nc;
Expand Down Expand Up @@ -107,13 +109,15 @@ void Plus_U::mix_locale(const UnitCell& ucell,

for (int T = 0; T < ucell.ntype; T++)
{
int target_l = get_orbital_corr(T);
if (target_l == -1)
continue;
if (!has_correlated_orbital(T))
{
continue;
}

for (int I = 0; I < ucell.atoms[T].na; I++)
{
const int iat = ucell.itia2iat(T, I);
int target_l = get_orbital_corr(T);

if (Plus_U::nspin == 4)
{
Expand Down
36 changes: 20 additions & 16 deletions source/source_lcao/module_dftu/dftu_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@ void Plus_U::cal_occ_pw(const int iter,
ModuleBase::timer::start("Plus_U", "cal_occ_pw");
this->copy_locale(cell);
this->zero_locale(cell);
const int nspin = PARAM.inp.nspin;
const int kpar = PARAM.inp.kpar;

if(this->device == "cpu")
{
Expand All @@ -37,7 +39,7 @@ void Plus_U::cal_occ_pw(const int iter,
const int npol = psi_p->get_npol();
for(int ik = 0; ik < psi_p->get_nk(); ik++)
{
int is = (Plus_U::nspin == 2) ? isk[ik] : 0;
int is = (nspin == 2) ? isk[ik] : 0;
psi_p->fix_k(ik);
onsite_p->tabulate_atomic(ik);

Expand All @@ -59,7 +61,7 @@ void Plus_U::cal_occ_pw(const int iter,
const int m_begin = target_l * target_l;
const int tlp1 = 2 * target_l + 1;
const int tlp1_2 = tlp1 * tlp1;
if(Plus_U::nspin == 4)
if(nspin == 4)
{
for(int ib = 0;ib<nbands;ib++)
{
Expand Down Expand Up @@ -117,7 +119,8 @@ void Plus_U::cal_occ_pw(const int iter,
const int npol = psi_p->get_npol();
for(int ik = 0; ik < psi_p->get_nk(); ik++)
{
int is = (Plus_U::nspin == 2) ? isk[ik] : 0;
int is = (nspin == 2) ? isk[ik] : 0;
const_cast<psi::Psi<std::complex<double>, base_device::DEVICE_GPU>*>(psi_p)->load_k_to_gpu(ik);
psi_p->fix_k(ik);
onsite_p->tabulate_atomic(ik);

Expand All @@ -138,7 +141,7 @@ void Plus_U::cal_occ_pw(const int iter,
const int m_begin = target_l * target_l;
const int tlp1 = 2 * target_l + 1;
const int tlp1_2 = tlp1 * tlp1;
if(Plus_U::nspin == 4)
if(nspin == 4)
{
for(int ib = 0;ib<nbands;ib++)
{
Expand Down Expand Up @@ -199,23 +202,23 @@ void Plus_U::cal_occ_pw(const int iter,
}
const int size = (2 * target_l + 1) * (2 * target_l + 1);

if(Plus_U::nspin != 4)
if(nspin != 4)
{
Parallel_Reduce::reduce_double_allpool(this->kpar,
Parallel_Reduce::reduce_double_allpool(kpar,
GlobalV::NPROC_IN_POOL,
this->locale[iat][target_l][0][0].c,
size);
if(Plus_U::nspin == 2)
if(nspin == 2)
{
Parallel_Reduce::reduce_double_allpool(this->kpar,
Parallel_Reduce::reduce_double_allpool(kpar,
GlobalV::NPROC_IN_POOL,
this->locale[iat][target_l][0][1].c,
size);
}
}
else
{
Parallel_Reduce::reduce_double_allpool(this->kpar,
Parallel_Reduce::reduce_double_allpool(kpar,
GlobalV::NPROC_IN_POOL,
this->locale[iat][target_l][0][0].c,
size * 4);
Expand All @@ -228,7 +231,7 @@ void Plus_U::cal_occ_pw(const int iter,
{
this->uom_array[eff_pot_pw_index[iat]+mm] = this->locale[iat][target_l][0][0].c[mm];
}
if(Plus_U::nspin == 2)
if(nspin == 2)
{
const int half_size = this->uom_array.size() / 2;
for(int mm=0;mm<size;mm++)
Expand All @@ -247,8 +250,8 @@ void Plus_U::cal_occ_pw(const int iter,
}

Plus_U::energy_u = 0.0;
const double weight_eu = (Plus_U::nspin == 1) ? 0.25 : (Plus_U::nspin == 2) ? 0.5 : 0.25;
const double diag_coeff = (Plus_U::nspin == 4) ? 1.0 : 0.5;
const double weight_eu = (nspin == 1) ? 1.0 : (nspin == 2) ? 0.5 : 0.25;
const double diag_coeff = (nspin == 4) ? 1.0 : 0.5;
// calculate VU and energy (locale already reduced above)
for(int iat = 0; iat < cell.nat; iat++)
{
Expand All @@ -265,7 +268,7 @@ void Plus_U::cal_occ_pw(const int iter,
std::complex<double>* vu_iat = &(this->eff_pot_pw[this->eff_pot_pw_index[iat]]);
const int m_size = 2 * target_l + 1;

if(Plus_U::nspin == 4)
if(nspin == 4)
{
for (int m1 = 0; m1 < m_size; m1++)
{
Expand Down Expand Up @@ -309,8 +312,8 @@ void Plus_U::cal_occ_pw(const int iter,
}
vu_iat[index[0]] = 0.5 * (vu_tmp[0] + vu_tmp[3]);
vu_iat[index[3]] = 0.5 * (vu_tmp[0] - vu_tmp[3]);
vu_iat[index[1]] = 0.5 * (vu_tmp[1] - std::complex<double>(0.0, 1.0) * vu_tmp[2]);
vu_iat[index[2]] = 0.5 * (vu_tmp[1] + std::complex<double>(0.0, 1.0) * vu_tmp[2]);
vu_iat[index[1]] = 0.5 * (vu_tmp[1] + std::complex<double>(0.0, 1.0) * vu_tmp[2]);
vu_iat[index[2]] = 0.5 * (vu_tmp[1] - std::complex<double>(0.0, 1.0) * vu_tmp[2]);
}
}
}
Expand All @@ -328,7 +331,7 @@ void Plus_U::cal_occ_pw(const int iter,
}
}
// spin-down channel for nspin=2
if(Plus_U::nspin == 2)
if(nspin == 2)
{
std::complex<double>* vu_iat1 = &(this->eff_pot_pw[this->eff_pot_pw.size()/2 + this->eff_pot_pw_index[iat]]);
for (int m1 = 0; m1 < m_size; m1++)
Expand All @@ -348,3 +351,4 @@ void Plus_U::cal_occ_pw(const int iter,
ModuleBase::timer::end("Plus_U", "cal_occ_pw");
}


16 changes: 1 addition & 15 deletions source/source_lcao/module_dftu/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,19 +1,5 @@
abacus_disable_feature_definitions(__CUDA)

AddTest(
TARGET dftu_pw_test
LIBS base device parameter
LIBS ${math_libs} base device parameter
SOURCES dftu_pw_test.cpp
)

AddTest(
TARGET dftu_core_test
LIBS base device
SOURCES dftu_core_test.cpp
)

AddTest(
TARGET dftu_operator_test
LIBS base device
SOURCES dftu_operator_test.cpp
)
Loading
Loading