Skip to content
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
3 changes: 0 additions & 3 deletions source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,6 @@ class PW_BASIS_K_GPU_TEST : public ::testing::Test
int distribution_type = 1;
bool xprime = false;
const int nks = 1;
ModuleBase::Vector3<double>* kvec_d;
kvec_d = new ModuleBase::Vector3<double>[nks];
kvec_d[0].set(0, 0, 0);
// init
const int mypool = 0;
const int key = 1;
Expand Down
4 changes: 4 additions & 0 deletions source/source_cell/module_symmetry/symm_rho.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,10 @@ void Symmetry::rho_symmetry( double *rho,
{
ModuleBase::timer::tick("Symmetry","rho_symmetry");

assert(nr1>0);
assert(nr2>0);
assert(nr3>0);

// allocate flag for each FFT grid.
bool* symflag = new bool[nr1 * nr2 * nr3];
for (int i=0; i<nr1*nr2*nr3; i++)
Expand Down
2 changes: 1 addition & 1 deletion source/source_esolver/esolver_of_tool.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,7 @@ double ESolver_OF::cal_mu(double* pphi, double* pdEdphi, double nelec)
* @brief Rotate and renormalize the direction |d>,
* make it orthogonal to phi (<d|phi> = 0), and <d|d> = nelec
*/
void ESolver_OF::adjust_direction()
void ESolver_OF::adjust_direction(void)
{
// filter the high frequency term in direction if of_full_pw = false
if (!PARAM.inp.of_full_pw)
Expand Down
90 changes: 62 additions & 28 deletions source/source_estate/module_charge/charge_mixing_preconditioner.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,21 +6,32 @@

void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
{
if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) {
return;
}
ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_recip");

if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1)
{
return;
}

ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip");

const int nspin = PARAM.inp.nspin;

double fac = 0.0;
double gg0 = 0.0;
double amin = 0.0;

/// consider a resize for mixing_angle
int resize_tmp = 1;
if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2;
}
int resize_tmp = 1;
if (nspin == 4 && this->mixing_angle > 0)
{
resize_tmp = 2;
}

/// implement Kerker for density and magnetization separately
for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
for (int is = 0; is < nspin / resize_tmp; ++is)
{
const int is_idx = is * this->rhopw->npw;
/// new mixing method only support nspin=2 not nspin=4
if (is >= 1)
{
Expand All @@ -29,10 +40,10 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
#ifdef __DEBUG
assert(is == 1); // make sure break works
#endif
double is_mag = PARAM.inp.nspin - 1;
double is_mag = nspin - 1;
//for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
//{
// drhog[is * this->rhopw->npw + ig] *= 1;
// drhog[is_idx + ig] *= 1;
//}
break;
}
Expand All @@ -46,32 +57,49 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
}

gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);

const double gg0_amin = this->mixing_gg0_min / amin;

#ifdef _OPENMP
#pragma omp parallel for schedule(static, 512)
#endif
for (int ig = 0; ig < this->rhopw->npw; ++ig)
{
double gg = this->rhopw->gg[ig];
double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
drhog[is * this->rhopw->npw + ig] *= filter_g;
double filter_g = std::max(gg / (gg + gg0), gg0_amin);
drhog[is_idx + ig] *= filter_g;
}
}

ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip");
return;
}

void Charge_Mixing::Kerker_screen_real(double* drhor)
{
if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) {
return;
}
/// consider a resize for mixing_angle
ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_real");

if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1)
{
return;
}

ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real");

const int nspin = PARAM.inp.nspin;
assert(nspin==1 || nspin==2 || nspin==4);

/// consider a resize for mixing_angle
int resize_tmp = 1;
if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2;
}
if (nspin == 4 && this->mixing_angle > 0)
{
resize_tmp = 2;
}

std::vector<std::complex<double>> drhog(this->rhopw->npw * PARAM.inp.nspin / resize_tmp);
std::vector<double> drhor_filter(this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp);
for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
std::vector<std::complex<double>> drhog(this->rhopw->npw * nspin / resize_tmp);
std::vector<double> drhor_filter(this->rhopw->nrxx * nspin / resize_tmp);

for (int is = 0; is < nspin / resize_tmp; ++is)
{
// Note after this process some G which is higher than Gmax will be filtered.
// Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it.
Expand All @@ -82,7 +110,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
double gg0 = 0.0;
double amin = 0.0;

for (int is = 0; is < PARAM.inp.nspin / resize_tmp; is++)
for (int is = 0; is < nspin / resize_tmp; is++)
{

if (is >= 1)
Expand All @@ -92,8 +120,8 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
#ifdef __DEBUG
assert(is == 1); /// make sure break works
#endif
double is_mag = PARAM.inp.nspin - 1;
if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { is_mag = 1;
double is_mag = nspin - 1;
if (nspin == 4 && this->mixing_angle > 0) { is_mag = 1;
}
for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
{
Expand All @@ -111,6 +139,9 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
}

gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);

const int is_idx = is * this->rhopw->npw;
const double gg0_amin = this->mixing_gg0_min / amin;
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 512)
#endif
Expand All @@ -120,24 +151,27 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
// I have not decided how to handle gg=0 part, will be changed in future
//if (gg == 0)
//{
// drhog[is * this->rhopw->npw + ig] *= 0;
// drhog[is_idx + ig] *= 0;
// continue;
//}
double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
drhog[is * this->rhopw->npw + ig] *= (1 - filter_g);
double filter_g = std::max(gg / (gg + gg0), gg0_amin);
drhog[is_idx + ig] *= (1 - filter_g);
}
}
/// inverse FT
for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
for (int is = 0; is < nspin / resize_tmp; ++is)
{
this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor_filter.data() + is * this->rhopw->nrxx);
}

#ifdef _OPENMP
#pragma omp parallel for schedule(static, 512)
#endif
for (int ir = 0; ir < this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp; ir++)
for (int ir = 0; ir < this->rhopw->nrxx * nspin / resize_tmp; ir++)
{
drhor[ir] -= drhor_filter[ir];
}

ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real");
return;
}
12 changes: 7 additions & 5 deletions source/source_estate/module_charge/charge_mixing_residual.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,13 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
{
ModuleBase::TITLE("Charge_Mixing", "get_drho");
ModuleBase::timer::tick("Charge_Mixing", "get_drho");
const int nspin = PARAM.inp.nspin;
assert(nspin==1 || nspin==2 || nspin==4);
double drho = 0.0;

if (PARAM.inp.scf_thr_type == 1)
{
for (int is = 0; is < PARAM.inp.nspin; ++is)
for (int is = 0; is < nspin; ++is)
{
ModuleBase::GlobalFunc::NOTE("Perform FFT on rho(r) to obtain rho(G).");
chr->rhopw->real2recip(chr->rho[is], chr->rhog[is]);
Expand All @@ -23,15 +25,15 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
}

ModuleBase::GlobalFunc::NOTE("Calculate the charge difference between rho(G) and rho_save(G)");
std::vector<std::complex<double>> drhog(PARAM.inp.nspin * this->rhopw->npw);
std::vector<std::complex<double>> drhog(nspin * this->rhopw->npw);
#ifdef _OPENMP
#pragma omp parallel for collapse(2) schedule(static, 512)
#endif
for (int is = 0; is < PARAM.inp.nspin; ++is)
for (int is = 0; is < nspin; ++is)
{
for (int ig = 0; ig < this->rhopw->npw; ig++)
{
drhog[is * rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig];
drhog[is * this->rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig];
}
}

Expand All @@ -42,7 +44,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
{
// Note: Maybe it is wrong.
// The inner_product_real function (L1-norm) is different from that (L2-norm) in mixing.
for (int is = 0; is < PARAM.inp.nspin; is++)
for (int is = 0; is < nspin; is++)
{
if (is != 0 && is != 3 && PARAM.globalv.domag_z)
{
Expand Down
Loading
Loading