Skip to content

Commit

Permalink
Refactor:Remove GlobalC::ucell in module_dftu (#5655)
Browse files Browse the repository at this point in the history 
* update ucell in force_gamma/force_k

* change GlobalC::ucell in the force_stress

* change GlobalC::ucell in the grid_init.cpp

* update ucell in hamilt_lcao

* update ucell in LCAO_alloacte

* update ucell in record_adj

* update ucell in spra_dh.cpp

* update ucell in wavefunc_in_pw

* fix bug in wavefunc

* change GlobalC::ucell in dftu.cpp

* update ucell in dftu_force

* update ucell in dftu_force

* change ucell in dftu_occup

* change ucell in dftu_yukawa

* change ucell in dftu_io

* [pre-commit.ci lite] apply automatic fixes

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>
  • Loading branch information
@A-006 @pre-commit-ci-lite
A-006 and pre-commit-ci-lite[bot] authored Dec 1, 2024
1 parent 5f1b128 commit 4a1750a
Show file tree
Hide file tree
Showing 10 changed files with 233 additions and 183 deletions.
8 changes: 4 additions & 4 deletions source/module_esolver/esolver_ks_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
GlobalC::dftu.set_dmr(dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
}
// Calculate U and J if Yukawa potential is used
GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx);
GlobalC::dftu.cal_slater_UJ(ucell,this->pelec->charge->rho, this->pw_rho->nrxx);
}

#ifdef __DEEPKS
Expand Down Expand Up @@ -875,11 +875,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
{
const std::vector<std::vector<TK>>& tmp_dm
= dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
ModuleDFTU::dftu_cal_occup_m(iter, tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
ModuleDFTU::dftu_cal_occup_m(iter, ucell,tmp_dm, this->kv, this->p_chgmix->get_mixing_beta(), this->p_hamilt);
}
GlobalC::dftu.cal_energy_correction(istep);
GlobalC::dftu.cal_energy_correction(ucell,istep);
}
GlobalC::dftu.output();
GlobalC::dftu.output(ucell);
}

// (7) for deepks, calculate delta_e
Expand Down
3 changes: 2 additions & 1 deletion source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -300,7 +300,8 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
}
if (PARAM.inp.dft_plus_u == 2)
{
GlobalC::dftu.force_stress(pelec,
GlobalC::dftu.force_stress(ucell,
pelec,
pv,
fsr, // mohan 2024-06-16
force_dftu,
Expand Down
9 changes: 5 additions & 4 deletions source/module_hamilt_lcao/hamilt_lcaodft/wavefunc_in_pw.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,8 +304,7 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const UnitCell& ucell,
/* for(int is_N = 0; is_N < 2; is_N++)*/ //for rotate base
for(int is_N = 0; is_N < 1; is_N++)
{
if(L==0 && is_N==1) { continue;
}
if(L==0 && is_N==1) { continue;}
if(ucell.atoms[it].ncpp.has_so)
{
const double j = std::abs(double(L+is_N) - 0.5);
Expand Down Expand Up @@ -370,9 +369,11 @@ void Wavefunc_in_pw::produce_local_basis_in_pw(const UnitCell& ucell,
for(int m = 0;m<2*L+1;m++)
{
const int lm = L*L +m;
if (iwall + 2 * L + 1 > ucell.natomwfc) {

if (iwall + 2 * L + 1 > ucell.natomwfc)
{
ModuleBase::WARNING_QUIT("this->wf.atomic_wfc()", "error: too many wfcs");
}
}
for (int ig = 0; ig < npw; ig++)
{
aux[ig] = sk[ig] * ylm(lm,ig) * chiaux[ig];
Expand Down
31 changes: 17 additions & 14 deletions source/module_hamilt_lcao/module_dftu/dftu.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,37 +179,38 @@ void DFTU::init(UnitCell& cell, // unitcell class
{
std::stringstream sst;
sst << "initial_onsite.dm";
this->read_occup_m(sst.str());
this->read_occup_m(cell,sst.str());
#ifdef __MPI
this->local_occup_bcast();
this->local_occup_bcast(cell);
#endif

initialed_locale = true;
this->copy_locale();
this->copy_locale(cell);
}
else
{
if (PARAM.inp.init_chg == "file")
{
std::stringstream sst;
sst << PARAM.globalv.global_out_dir << "onsite.dm";
this->read_occup_m(sst.str());
this->read_occup_m(cell,sst.str());
#ifdef __MPI
this->local_occup_bcast();
this->local_occup_bcast(cell);
#endif
initialed_locale = true;
}
else
{
this->zero_locale();
this->zero_locale(cell);
}
}

ModuleBase::Memory::record("DFTU::locale", sizeof(double) * num_locale);
return;
}

void DFTU::cal_energy_correction(const int istep)
void DFTU::cal_energy_correction(const UnitCell& ucell,
const int istep)
{
ModuleBase::TITLE("DFTU", "cal_energy_correction");
ModuleBase::timer::tick("DFTU", "cal_energy_correction");
Expand All @@ -221,18 +222,18 @@ void DFTU::cal_energy_correction(const int istep)
this->EU = 0.0;
double EU_dc = 0.0;

for (int T = 0; T < GlobalC::ucell.ntype; T++)
for (int T = 0; T < ucell.ntype; T++)
{
const int NL = GlobalC::ucell.atoms[T].nwl + 1;
const int NL = ucell.atoms[T].nwl + 1;
const int LC = orbital_corr[T];
for (int I = 0; I < GlobalC::ucell.atoms[T].na; I++)
for (int I = 0; I < ucell.atoms[T].na; I++)
{
if (LC == -1)
{
continue;
}

const int iat = GlobalC::ucell.itia2iat(T, I);
const int iat = ucell.itia2iat(T, I);
const int L = orbital_corr[T];

for (int l = 0; l < NL; l++)
Expand All @@ -242,7 +243,7 @@ void DFTU::cal_energy_correction(const int istep)
continue;
}

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

const int m_tot = 2 * l + 1;

Expand Down Expand Up @@ -422,22 +423,24 @@ const hamilt::HContainer<double>* DFTU::get_dmr(int ispin) const
//! dftu occupation matrix for gamma only using dm(double)
template <>
void dftu_cal_occup_m(const int iter,
const UnitCell& ucell,
const std::vector<std::vector<double>>& dm,
const K_Vectors& kv,
const double& mixing_beta,
hamilt::Hamilt<double>* p_ham)
{
GlobalC::dftu.cal_occup_m_gamma(iter, dm, mixing_beta, p_ham);
GlobalC::dftu.cal_occup_m_gamma(iter, ucell ,dm, mixing_beta, p_ham);
}

//! dftu occupation matrix for multiple k-points using dm(complex)
template <>
void dftu_cal_occup_m(const int iter,
const UnitCell& ucell,
const std::vector<std::vector<std::complex<double>>>& dm,
const K_Vectors& kv,
const double& mixing_beta,
hamilt::Hamilt<std::complex<double>>* p_ham)
{
GlobalC::dftu.cal_occup_m_k(iter, dm, kv, mixing_beta, p_ham);
GlobalC::dftu.cal_occup_m_k(iter,ucell, dm, kv, mixing_beta, p_ham);
}
} // namespace ModuleDFTU
90 changes: 53 additions & 37 deletions source/module_hamilt_lcao/module_dftu/dftu.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ class DFTU
);

// calculate the energy correction
void cal_energy_correction(const int istep);
void cal_energy_correction(const UnitCell& ucell, const int istep);
double get_energy(){return EU;}
void uramping_update(); // update U by uramping
bool u_converged(); // check if U is converged
Expand Down Expand Up @@ -89,16 +89,25 @@ class DFTU
//=============================================================
public:
// calculate the local occupation number matrix
void cal_occup_m_k(const int iter, const std::vector<std::vector<std::complex<double>>>& dm_k, const K_Vectors& kv, const double& mixing_beta, hamilt::Hamilt<std::complex<double>>* p_ham);
void cal_occup_m_gamma(const int iter, const std::vector<std::vector<double>>& dm_gamma, const double& mixing_beta, hamilt::Hamilt<double>* p_ham);
void cal_occup_m_k(const int iter,
const UnitCell& ucell,
const std::vector<std::vector<std::complex<double>>>& dm_k,
const K_Vectors& kv,
const double& mixing_beta,
hamilt::Hamilt<std::complex<double>>* p_ham);
void cal_occup_m_gamma(const int iter,
const UnitCell& ucell,
const std::vector<std::vector<double>>& dm_gamma,
const double& mixing_beta,
hamilt::Hamilt<double>* p_ham);

// dftu can be calculated only after locale has been initialed
bool initialed_locale = false;

private:
void copy_locale();
void zero_locale();
void mix_locale(const double& mixing_beta);
void copy_locale(const UnitCell& ucell);
void zero_locale(const UnitCell& ucell);
void mix_locale(const UnitCell& ucell,const double& mixing_beta);

public:
// local occupancy matrix of the correlated subspace
Expand Down Expand Up @@ -147,13 +156,14 @@ class DFTU
// dim = 4-6 : dS * dR, for stress

void folding_matrix_k(
const UnitCell &ucell,
ForceStressArrays &fsr,
const Parallel_Orbitals &pv,
const int ik,
const int dim1,
const int dim2,
std::complex<double>* mat_k,
const std::vector<ModuleBase::Vector3<double>> &kvec_d);
const int ik,
const int dim1,
const int dim2,
std::complex<double>* mat_k,
const std::vector<ModuleBase::Vector3<double>> &kvec_d);


/**
Expand All @@ -169,38 +179,40 @@ class DFTU
//=============================================================
public:

void force_stress(const elecstate::ElecState* pelec,
const Parallel_Orbitals& pv,
ForceStressArrays& fsr,
ModuleBase::matrix& force_dftu,
ModuleBase::matrix& stress_dftu,
const K_Vectors& kv);
void force_stress(const UnitCell& ucell,
const elecstate::ElecState* pelec,
const Parallel_Orbitals& pv,
ForceStressArrays& fsr,
ModuleBase::matrix& force_dftu,
ModuleBase::matrix& stress_dftu,
const K_Vectors& kv);

private:

void cal_force_k(
ForceStressArrays &fsr,
const Parallel_Orbitals &pv,
const int ik,
const std::complex<double>* rho_VU,
ModuleBase::matrix& force_dftu,
const std::vector<ModuleBase::Vector3<double>>& kvec_d);
void cal_force_k(const UnitCell &ucell,
ForceStressArrays &fsr,
const Parallel_Orbitals &pv,
const int ik,
const std::complex<double>* rho_VU,
ModuleBase::matrix& force_dftu,
const std::vector<ModuleBase::Vector3<double>>& kvec_d);

void cal_stress_k(
const UnitCell &ucell,
ForceStressArrays &fsr,
const Parallel_Orbitals &pv,
const int ik,
const std::complex<double>* rho_VU,
ModuleBase::matrix& stress_dftu,
const std::vector<ModuleBase::Vector3<double>>& kvec_d);

void cal_force_gamma(
const double* rho_VU,
const Parallel_Orbitals &pv,
double* dsloc_x,
double* dsloc_y,
double* dsloc_z,
ModuleBase::matrix& force_dftu);
void cal_force_gamma(const UnitCell &ucell,
const double* rho_VU,
const Parallel_Orbitals &pv,
double* dsloc_x,
double* dsloc_y,
double* dsloc_z,
ModuleBase::matrix& force_dftu);

void cal_stress_gamma(
const UnitCell &ucell,
Expand All @@ -218,12 +230,15 @@ class DFTU
// For reading/writing/broadcasting/copying relevant data structures
//=============================================================
public:
void output();
void output(const UnitCell& ucell);

private:
void write_occup_m(std::ofstream& ofs, bool diag=false);
void read_occup_m(const std::string& fn);
void local_occup_bcast();
void write_occup_m(const UnitCell& ucell,
std::ofstream& ofs,
bool diag=false);
void read_occup_m(const UnitCell& ucell,
const std::string& fn);
void local_occup_bcast(const UnitCell& ucell);

//=============================================================
// In dftu_yukawa.cpp
Expand All @@ -232,15 +247,15 @@ class DFTU

public:
bool Yukawa; // 1:use Yukawa potential; 0: do not use Yukawa potential
void cal_slater_UJ(double** rho, const int& nrxx);
void cal_slater_UJ(const UnitCell& ucell, double** rho, const int& nrxx);

private:
double lambda; // the parameter in Yukawa potential
std::vector<std::vector<std::vector<std::vector<double>>>> Fk; // slater integral:Fk[T][L][N][k]
std::vector<std::vector<std::vector<double>>> U_Yukawa; // U_Yukawa[T][L][N]
std::vector<std::vector<std::vector<double>>> J_Yukawa; // J_Yukawa[T][L][N]

void cal_slater_Fk(const int L, const int T); // L:angular momnet, T:atom type
void cal_slater_Fk(const UnitCell& ucell,const int L, const int T); // L:angular momnet, T:atom type
void cal_yukawa_lambda(double** rho, const int& nrxx);

double spherical_Bessel(const int k, const double r, const double lambda);
Expand All @@ -267,6 +282,7 @@ class DFTU

template <typename T>
void dftu_cal_occup_m(const int iter,
const UnitCell& ucell,
const std::vector<std::vector<T>>& dm,
const K_Vectors& kv,
const double& mixing_beta,
Expand Down
25 changes: 13 additions & 12 deletions source/module_hamilt_lcao/module_dftu/dftu_folding.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,6 +126,7 @@ void DFTU::fold_dSR_gamma(
}

void DFTU::folding_matrix_k(
const UnitCell &ucell,
ForceStressArrays &fsr,
const Parallel_Orbitals &pv,
const int ik,
Expand Down Expand Up @@ -153,26 +154,26 @@ void DFTU::folding_matrix_k(
ModuleBase::Vector3<double> dtau2;
ModuleBase::Vector3<double> tau0;

for (int T1 = 0; T1 < GlobalC::ucell.ntype; ++T1)
for (int T1 = 0; T1 < ucell.ntype; ++T1)
{
Atom* atom1 = &GlobalC::ucell.atoms[T1];
Atom* atom1 = &ucell.atoms[T1];
for (int I1 = 0; I1 < atom1->na; ++I1)
{
tau1 = atom1->tau[I1];
GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1);
Atom* atom1 = &GlobalC::ucell.atoms[T1];
const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0);
GlobalC::GridD.Find_atom(ucell, tau1, T1, I1);
Atom* atom1 = &ucell.atoms[T1];
const int start1 = ucell.itiaiw2iwt(T1, I1, 0);

// (2) search among all adjacent atoms.
for (int ad = 0; ad < GlobalC::GridD.getAdjacentNum() + 1; ++ad)
{
const int T2 = GlobalC::GridD.getType(ad);
const int I2 = GlobalC::GridD.getNatom(ad);
Atom* atom2 = &GlobalC::ucell.atoms[T2];
Atom* atom2 = &ucell.atoms[T2];

tau2 = GlobalC::GridD.getAdjacentTau(ad);
dtau = tau2 - tau1;
double distance = dtau.norm() * GlobalC::ucell.lat0;
double distance = dtau.norm() * ucell.lat0;
double rcut = orb_cutoff_[T1] + orb_cutoff_[T2];

bool adj = false;
Expand All @@ -192,11 +193,11 @@ void DFTU::folding_matrix_k(
dtau1 = tau0 - tau1;
dtau2 = tau0 - tau2;

double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
double distance1 = dtau1.norm() * ucell.lat0;
double distance2 = dtau2.norm() * ucell.lat0;

double rcut1 = orb_cutoff_[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
double rcut2 = orb_cutoff_[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
double rcut1 = orb_cutoff_[T1] + ucell.infoNL.Beta[T0].get_rcut_max();
double rcut2 = orb_cutoff_[T2] + ucell.infoNL.Beta[T0].get_rcut_max();

if (distance1 < rcut1 && distance2 < rcut2)
{
Expand All @@ -209,7 +210,7 @@ void DFTU::folding_matrix_k(
if (adj)
{
// (3) calculate the nu of atom (T2, I2)
const int start2 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0);
const int start2 = ucell.itiaiw2iwt(T2, I2, 0);
//------------------------------------------------
// exp(k dot dR)
// dR is the index of box in Crystal coordinates
Expand Down
Loading

0 comments on commit 4a1750a

Please sign in to comment.