Skip to content

add <psi|r|beta> function #8

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
ESROAMER merged 1 commit into from
Mar 12, 2025
Merged

add <psi|r|beta> function #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
3 changes: 3 additions & 0 deletions source/module_basis/module_ao/ORB_nonlocal_lm.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,9 @@ class Numerical_Nonlocal_Lm
const double& getKpoint(const int &ik) const { return this->k_radial[ik]; }
const double* getBeta_k() const { return this->beta_k; }
const double& getBeta_k(const int &ik) const { return this->beta_k[ik]; }

const int& getNk() const { return nk; }
const double& getDruniform() const { return dr_uniform; }

// enables deep copy
Numerical_Nonlocal_Lm& operator= (const Numerical_Nonlocal_Lm& nol );
Expand Down
306 changes: 306 additions & 0 deletions source/module_io/cal_r_overlap_R.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "module_base/timer.h"
#include "module_cell/module_neighbor/sltk_grid_driver.h"
#include "module_hamilt_pw/hamilt_pwdft/global.h"
#include "module_base/mathzone_add1.h"

cal_r_overlap_R::cal_r_overlap_R()
{
Expand Down Expand Up @@ -224,6 +225,230 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const LCAO_Orbitals& orb)
}
}

void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const LCAO_Orbitals& orb)
{
const InfoNonlocal& infoNL_ = GlobalC::ucell.infoNL;

int orb_r_ntype = 0;
int mat_Nr = orb.Phi[0].PhiLN(0, 0).getNr();
int count_Nr = 0;

orbs.resize(orb.get_ntype());
for (int T = 0; T < orb.get_ntype(); ++T)
{
count_Nr = orb.Phi[T].PhiLN(0, 0).getNr();
if (count_Nr > mat_Nr)
{
mat_Nr = count_Nr;
orb_r_ntype = T;
}

orbs[T].resize(orb.Phi[T].getLmax() + 1);
for (int L = 0; L <= orb.Phi[T].getLmax(); ++L)
{
orbs[T][L].resize(orb.Phi[T].getNchi(L));
for (int N = 0; N < orb.Phi[T].getNchi(L); ++N)
{
const auto& orb_origin = orb.Phi[T].PhiLN(L, N);
orbs[T][L][N].set_orbital_info(orb_origin.getLabel(),
orb_origin.getType(),
orb_origin.getL(),
orb_origin.getChi(),
orb_origin.getNr(),
orb_origin.getRab(),
orb_origin.getRadial(),
Numerical_Orbital_Lm::Psi_Type::Psi,
orb_origin.getPsi(),
static_cast<int>(orb_origin.getNk() * kmesh_times) | 1,
orb_origin.getDk(),
orb_origin.getDruniform(),
false,
true,
PARAM.inp.cal_force);
}
}
}

orb_r.set_orbital_info(orbs[orb_r_ntype][0][0].getLabel(), // atom label
orb_r_ntype, // atom type
1, // angular momentum L
1, // number of orbitals of this L , just N
orbs[orb_r_ntype][0][0].getNr(), // number of radial mesh
orbs[orb_r_ntype][0][0].getRab(), // the mesh interval in radial mesh
orbs[orb_r_ntype][0][0].getRadial(), // radial mesh value(a.u.)
Numerical_Orbital_Lm::Psi_Type::Psi,
orbs[orb_r_ntype][0][0].getRadial(), // radial wave function
orbs[orb_r_ntype][0][0].getNk(),
orbs[orb_r_ntype][0][0].getDk(),
orbs[orb_r_ntype][0][0].getDruniform(),
false,
true,
PARAM.inp.cal_force);

orbs_nonlocal.resize(orb.get_ntype());
for (int T = 0; T < orb.get_ntype(); ++T)
{
const int nproj = infoNL_.nproj[T];
orbs_nonlocal[T].resize(nproj);
for (int ip = 0; ip < nproj; ip++)
{
int nr = infoNL_.Beta[T].Proj[ip].getNr();
double dr_uniform = 0.01;
int nr_uniform = static_cast<int>((infoNL_.Beta[T].Proj[ip].getRadial(nr-1) - infoNL_.Beta[T].Proj[ip].getRadial(0))/dr_uniform) + 1;
double* rad = new double[nr_uniform];
double* rab = new double[nr_uniform];
for (int ir = 0; ir < nr_uniform; ir++)
{
rad[ir] = ir*dr_uniform;
rab[ir] = dr_uniform;
}
double* y2 = new double[nr];
double* Beta_r_uniform = new double[nr_uniform];
double* dbeta_uniform = new double[nr_uniform];
ModuleBase::Mathzone_Add1::SplineD2(infoNL_.Beta[T].Proj[ip].getRadial(), infoNL_.Beta[T].Proj[ip].getBeta_r(), nr, 0.0, 0.0, y2);
ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(
infoNL_.Beta[T].Proj[ip].getRadial(),
infoNL_.Beta[T].Proj[ip].getBeta_r(),
y2,
nr,
rad,
nr_uniform,
Beta_r_uniform,
dbeta_uniform
);

// linear extrapolation at the zero point
if (infoNL_.Beta[T].Proj[ip].getRadial(0) > 1e-10)
{
double slope = (infoNL_.Beta[T].Proj[ip].getBeta_r(1) - infoNL_.Beta[T].Proj[ip].getBeta_r(0)) / (infoNL_.Beta[T].Proj[ip].getRadial(1) - infoNL_.Beta[T].Proj[ip].getRadial(0));
Beta_r_uniform[0] = infoNL_.Beta[T].Proj[ip].getBeta_r(0) - slope * infoNL_.Beta[T].Proj[ip].getRadial(0);
}

orbs_nonlocal[T][ip].set_orbital_info(infoNL_.Beta[T].getLabel(),
infoNL_.Beta[T].getType(),
infoNL_.Beta[T].Proj[ip].getL(),
1,
nr_uniform,
rab,
rad,
Numerical_Orbital_Lm::Psi_Type::Psi,
Beta_r_uniform,
static_cast<int>(infoNL_.Beta[T].Proj[ip].getNk() * kmesh_times) | 1,
infoNL_.Beta[T].Proj[ip].getDk(),
infoNL_.Beta[T].Proj[ip].getDruniform(),
false,
true,
PARAM.inp.cal_force);

delete [] rad;
delete [] rab;
delete [] y2;
delete [] Beta_r_uniform;
delete [] dbeta_uniform;
}
}

for (int TA = 0; TA < orb.get_ntype(); ++TA)
{
for (int TB = 0; TB < orb.get_ntype(); ++TB)
{
for (int LA = 0; LA <= orb.Phi[TA].getLmax(); ++LA)
{
for (int NA = 0; NA < orb.Phi[TA].getNchi(LA); ++NA)
{
for (int ip = 0; ip < infoNL_.nproj[TB]; ip++)
{
center2_orb11_nonlocal[TA][TB][LA][NA].insert(
std::make_pair(ip, Center2_Orb::Orb11(orbs[TA][LA][NA], orbs_nonlocal[TB][ip], psb_, MGT)));
}
}
}
}
}

for (int TA = 0; TA < orb.get_ntype(); ++TA)
{
for (int TB = 0; TB < orb.get_ntype(); ++TB)
{
for (int LA = 0; LA <= orb.Phi[TA].getLmax(); ++LA)
{
for (int NA = 0; NA < orb.Phi[TA].getNchi(LA); ++NA)
{
for (int ip = 0; ip < infoNL_.nproj[TB]; ip++)
{
center2_orb21_r_nonlocal[TA][TB][LA][NA].insert(std::make_pair(
ip,
Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs_nonlocal[TB][ip], psb_, MGT)));
}
}
}
}
}

for (auto& co1: center2_orb11_nonlocal)
{
for (auto& co2: co1.second)
{
for (auto& co3: co2.second)
{
for (auto& co4: co3.second)
{
for (auto& co5: co4.second)
{
co5.second.init_radial_table();
}
}
}
}
}

for (auto& co1: center2_orb21_r_nonlocal)
{
for (auto& co2: co1.second)
{
for (auto& co3: co2.second)
{
for (auto& co4: co3.second)
{
for (auto& co5: co4.second)
{
co5.second.init_radial_table();
}
}
}
}
}

iw2it.resize(PARAM.globalv.nlocal);
iw2ia.resize(PARAM.globalv.nlocal);
iw2iL.resize(PARAM.globalv.nlocal);
iw2iN.resize(PARAM.globalv.nlocal);
iw2im.resize(PARAM.globalv.nlocal);

int iw = 0;
for (int it = 0; it < GlobalC::ucell.ntype; it++)
{
for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++)
{
for (int iL = 0; iL < GlobalC::ucell.atoms[it].nwl + 1; iL++)
{
for (int iN = 0; iN < GlobalC::ucell.atoms[it].l_nchi[iL]; iN++)
{
for (int im = 0; im < (2 * iL + 1); im++)
{
iw2it[iw] = it;
iw2ia[iw] = ia;
iw2iL[iw] = iL;
iw2iN[iw] = iN;
iw2im[iw] = im;
iw++;
}
}
}
}
}
}

void cal_r_overlap_R::init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb)
{
ModuleBase::TITLE("cal_r_overlap_R", "init");
Expand All @@ -236,6 +461,21 @@ void cal_r_overlap_R::init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb
ModuleBase::timer::tick("cal_r_overlap_R", "init");
return;
}

void cal_r_overlap_R::init_nonlocal(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb)
{
ModuleBase::TITLE("cal_r_overlap_R", "init_nonlocal");
ModuleBase::timer::tick("cal_r_overlap_R", "init_nonlocal");
this->ParaV = &pv;

initialize_orb_table(orb);
construct_orbs_and_nonlocal_and_orb_r(orb);

ModuleBase::timer::tick("cal_r_overlap_R", "init_nonlocal");
return;
}


ModuleBase::Vector3<double> cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
Expand Down Expand Up @@ -279,6 +519,72 @@ ModuleBase::Vector3<double> cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vec

return temp_prp;
}

void cal_r_overlap_R::get_psi_r_beta(std::vector<ModuleBase::Vector3<double>>& nlm,
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2)
{
ModuleBase::Vector3<double> origin_point(0.0, 0.0, 0.0);
double factor = sqrt(ModuleBase::FOUR_PI / 3.0);
const ModuleBase::Vector3<double>& distance = R2 - R1;
const InfoNonlocal& infoNL_ = GlobalC::ucell.infoNL;
const int nproj = infoNL_.nproj[T2];
if (nproj == 0)
{
nlm.resize(1);
return;
}

int natomwfc = 0;
for (int ip = 0; ip < nproj; ip++)
{
const int L2 = infoNL_.Beta[T2].Proj[ip].getL(); // mohan add 2021-05-07
natomwfc += 2 * L2 + 1;
}
nlm.resize(natomwfc);

int index = 0;
for (int ip = 0; ip < nproj; ip++)
{
const int L2 = infoNL_.Beta[T2].Proj[ip].getL();
for (int m2 = 0; m2 < 2 * L2 + 1; m2++)
{
double overlap_o
= center2_orb11_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, distance, m1, m2);

double overlap_x = -1 * factor
* center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point,
distance,
m1,
1,
m2); // m = 1

double overlap_y = -1 * factor
* center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point,
distance,
m1,
2,
m2); // m = -1

double overlap_z = factor
* center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point,
distance,
m1,
0,
m2); // m = 0

nlm[index] = ModuleBase::Vector3<double>(overlap_x, overlap_y, overlap_z) + R1 * overlap_o;
index++;
}
}
}


void cal_r_overlap_R::out_rR(const int& istep)
{
ModuleBase::TITLE("cal_r_overlap_R", "out_rR");
Expand Down
14 changes: 13 additions & 1 deletion source/module_io/cal_r_overlap_R.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ class cal_r_overlap_R
bool binary = false;

void init(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
ModuleBase::Vector3<double> get_psi_r_psi(
void init_nonlocal(const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
ModuleBase::Vector3<double> get_psi_r_psi(
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
Expand All @@ -43,13 +44,24 @@ class cal_r_overlap_R
const int& m2,
const int& N2
);
void get_psi_r_beta(
std::vector<ModuleBase::Vector3<double>>& nlm,
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2
);
void out_rR(const int& istep);
void out_rR_other(const int& istep, const std::set<Abfs::Vector3_Order<int>>& output_R_coor);


private:
void initialize_orb_table(const LCAO_Orbitals& orb);
void construct_orbs_and_orb_r(const LCAO_Orbitals& orb);
void construct_orbs_and_nonlocal_and_orb_r(const LCAO_Orbitals& orb);

std::vector<int> iw2ia;
std::vector<int> iw2iL;
Expand Down