Skip to content

Fix: set the maximal cutoff radius of orbital as the lower bound of search_radius in neighboring search tasks carried in out_mat_l #6653

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

Merged
merged 3 commits into from
Oct 29, 2025
Merged

Fix: set the maximal cutoff radius of orbital as the lower bound of search_radius in neighboring search tasks carried in out_mat_l #6653

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
47975c1
Fix: adjust neighbor search loop condition in AngularMomentumCalculator
kirk0830 Oct 29, 2025
59ead16
Fix: adjust search_radius handling and improve neighbor search in Ang…
kirk0830 Oct 29, 2025
4d5e538
Merge branch 'develop' into fix-nadj-calcpLp
kirk0830 Oct 29, 2025
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
50 changes: 32 additions & 18 deletions source/source_io/cal_pLpR.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,12 +150,21 @@ ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator(

// for neighbor list search
double temp = -1.0;
if (search_radius < rcut_max)
{
*ofs_ << "Find the `search_radius` from the input file being smaller than the \n"
"`rcut_max` of the orbitals.\n"
<< "Reset the `search_radius` (" << search_radius << ") "
<< "to `rcut_max` ("<< rcut_max << ")."
<< std::endl;
// we don't really set, but use std::max to mask :)
}
temp = atom_arrange::set_sr_NL(*ofs_,
PARAM.inp.out_level,
search_radius,
std::max(search_radius, rcut_max),
ucell.infoNL.get_rcutmax_Beta(),
PARAM.globalv.gamma_only_local);
temp = std::max(temp, search_radius);
temp = std::max(temp, std::max(search_radius, rcut_max));
this->neighbor_searcher_ = std::unique_ptr<Grid_Driver>(new Grid_Driver(tdestructor, tgrid));
atom_arrange::search(searchpbc,
*ofs_,
Expand Down Expand Up @@ -190,24 +199,29 @@ void ModuleIO::AngularMomentumCalculator::kernel(
fmtstr += "%" + std::to_string(precision*2) + "." + std::to_string(precision) + "e\n";
FmtCore fmt(fmtstr);

ModuleBase::Vector3<double> ri, rj, dr;
// placeholders
std::complex<double> val = 0;
ModuleBase::Vector3<double> taui; // the origin position
ModuleBase::Vector3<double> dtau; // the displacement
AdjacentAtomInfo adjinfo; // adjacent atom information carrier
for (int it = 0; it < ucell.ntype; it++)
{
const Atom& atyp_i = ucell.atoms[it];
for (int ia = 0; ia < atyp_i.na; ia++)
{
ri = atyp_i.tau[ia];
neighbor_searcher_->Find_atom(ucell, ri, it, ia);
for (int ia_adj = 0; ia_adj < neighbor_searcher_->getAdjacentNum(); ia_adj++)
taui = ucell.get_tau(ucell.itia2iat(it, ia));
neighbor_searcher_->Find_atom(ucell, taui, it, ia, &adjinfo);
for (int ia_adj = 0; ia_adj < adjinfo.adj_num + 1; ia_adj++) // "+1" is to include itself
{
rj = neighbor_searcher_->getAdjacentTau(ia_adj);
int jt = neighbor_searcher_->getType(ia_adj);
int jt = adjinfo.ntype[ia_adj]; // ityp
int ja = adjinfo.natom[ia_adj]; // iat with in atomtype
const Atom& atyp_j = ucell.atoms[jt];
int ja = neighbor_searcher_->getNatom(ia_adj);
dr = (ri - rj) * ucell.lat0;
const ModuleBase::Vector3<int> iR = neighbor_searcher_->getBox(ia_adj);
// the two-center-integral
const ModuleBase::Vector3<int> iR = adjinfo.box[ia_adj];
dtau = ucell.cal_dtau(ucell.itia2iat(it, ia),
ucell.itia2iat(jt, ja),
iR) * ucell.lat0; // convert to unit of Bohr

// nested loop: calculate the two-center-integral
for (int li = 0; li < atyp_i.nwl + 1; li++)
{
for (int iz = 0; iz < atyp_i.l_nchi[li]; iz++)
Expand All @@ -220,21 +234,20 @@ void ModuleIO::AngularMomentumCalculator::kernel(
{
for (int mj = -lj; mj <= lj; mj++)
{
std::complex<double> val = 0;
if (dir == 'x')
{
val = cal_LxijR(calculator_,
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
}
else if (dir == 'y')
{
val = cal_LyijR(calculator_,
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
}
else if (dir == 'z')
{
val = cal_LzijR(calculator_,
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dr);
it, ia, li, iz, mi, jt, ja, lj, jz, mj, dtau);
}

*ofs << fmt.format(
Expand Down Expand Up @@ -266,7 +279,7 @@ void ModuleIO::AngularMomentumCalculator::calculate(
}
std::ofstream ofout;
const std::string dir = "xyz";
const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm <a|L|b>\n"
const std::string title = "# it ia il iz im iRx iRy iRz jt ja jl jz jm Re[<a|L|b>] Im[<a|L|b>]\n"
"# it: atomtype index of the first atom\n"
"# ia: atomic index of the first atom within the atomtype\n"
"# il: angular momentum index of the first atom\n"
Expand All @@ -278,7 +291,8 @@ void ModuleIO::AngularMomentumCalculator::calculate(
"# jl: angular momentum index of the second atom\n"
"# jz: zeta function index of the second atom\n"
"# jm: magnetic quantum number of the second atom\n"
"# <a|L|b>: the value of the matrix element\n";
"# Re[<a|L|b>], Im[<a|L|b>]: the real and imaginary parts "
"of the value of the matrix element\n";

for (char d : dir)
{
Expand Down
Loading