Skip to content

Feature: mixing density and magnetism {mx,my,mz} separately for non-colinear calculation #3330

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 8 commits into from
Dec 14, 2023
Merged

Feature: mixing density and magnetism {mx,my,mz} separately for non-colinear calculation #3330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e6a7403
new mixing for nspin=4
WHUweiqingzhou Dec 11, 2023
6e77975
revert get_drho
WHUweiqingzhou Dec 12, 2023
3c27580
change default value and open mixing_beta_mag for nspin=4
WHUweiqingzhou Dec 13, 2023
dc300bc
support new mixing of nspin=4 also for PW basis
WHUweiqingzhou Dec 13, 2023
b57fae7
remove seting nspin in front of seting mixing in Default_2()
WHUweiqingzhou Dec 13, 2023
78a790d
update ref of 104_PW_NC_magnetic
WHUweiqingzhou Dec 13, 2023
05a97a3
update INPUT of 204_NO_KP_NC
WHUweiqingzhou Dec 13, 2023
3d04a10
Merge branch 'develop' into mixing_nspin4
WHUweiqingzhou Dec 13, 2023
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
113 changes: 90 additions & 23 deletions source/module_elecstate/module_charge/charge_mixing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,9 +38,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
GlobalV::ofs_running<<"mixing_beta: "<< this->mixing_beta <<std::endl;
GlobalV::ofs_running<<"mixing_gg0: "<< this->mixing_gg0 <<std::endl;
GlobalV::ofs_running<<"mixing_gg0_min: "<< GlobalV::MIXING_GG0_MIN <<std::endl;
if (GlobalV::NSPIN==2)
if (GlobalV::NSPIN==2 || GlobalV::NSPIN==4)
{
GlobalV::ofs_running<<"mixing_beta_mag: "<< GlobalV::MIXING_BETA_MAG <<std::endl;
GlobalV::ofs_running<<"mixing_beta_mag: "<< this->mixing_beta_mag <<std::endl;
GlobalV::ofs_running<<"mixing_gg0_mag: "<< GlobalV::MIXING_GG0_MAG <<std::endl;
}
GlobalV::ofs_running<<"mixing_ndim: "<< this->mixing_ndim <<std::endl;
Expand Down Expand Up @@ -271,7 +271,6 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)

void Charge_Mixing::mix_rho_recip_new(Charge* chr)
{
// not support nspin=4 yet 2023/11/17
// old support see mix_rho_recip()
if (GlobalV::double_grid)
{
Expand All @@ -281,7 +280,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
std::complex<double>* rhog_in = nullptr;
std::complex<double>* rhog_out = nullptr;

if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
if (GlobalV::NSPIN == 1)
{
rhog_in = chr->rhog_save[0];
rhog_out = chr->rhog[0];
Expand Down Expand Up @@ -316,8 +315,38 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
};
this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true);
}
else if (GlobalV::NSPIN == 4)
{
// I do not merge this part with nspin=2 because the method of nspin=2 is almost done,
// while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future.
rhog_in = chr->rhog_save[0];
rhog_out = chr->rhog[0];
const int npw = this->rhopw->npw;
auto screen = std::bind(&Charge_Mixing::Kerker_screen_recip_new, this, std::placeholders::_1); // use old one
auto twobeta_mix
= [this, npw](std::complex<double>* out, const std::complex<double>* in, const std::complex<double>* sres) {
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = 0; i < npw; ++i)
{
out[i] = in[i] + this->mixing_beta * sres[i];
}
// magnetism, mx, my, mz
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = npw; i < 4 * npw; ++i)
{
out[i] = in[i] + this->mixing_beta_mag * sres[i];
}
};
this->mixing->push_data(this->rho_mdata, rhog_in, rhog_out, screen, twobeta_mix, true);
}

// can choose inner_product_recip_new1 or inner_product_recip_new2
// inner_product_recip_new1 is a simple sum
// inner_product_recip_new2 is a hartree-like sum, unit is Ry
auto inner_product_new
= std::bind(&Charge_Mixing::inner_product_recip_new2, this, std::placeholders::_1, std::placeholders::_2);
auto inner_product_old
Expand All @@ -326,7 +355,7 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
{
this->mixing->cal_coef(this->rho_mdata, inner_product_new);
}
else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
else if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4) // nspin=4 can use old inner_product
{
this->mixing->cal_coef(this->rho_mdata, inner_product_old);
}
Expand All @@ -344,6 +373,8 @@ void Charge_Mixing::mix_rho_recip_new(Charge* chr)
{
chr->rhopw->recip2real(chr->rhog[is], chr->rho[is]);
}

// renormalize rho in R-space would induce a error in K-space
//chr->renormalize_rho();

// For kinetic energy density
Expand Down Expand Up @@ -400,7 +431,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
{
double* rhor_in;
double* rhor_out;
if (GlobalV::NSPIN == 1 || GlobalV::NSPIN == 4)
if (GlobalV::NSPIN == 1)
{
rhor_in = chr->rho_save[0];
rhor_out = chr->rho[0];
Expand All @@ -420,19 +451,47 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = 0; i < nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta * sres[i];
}
for (int i = 0; i < nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta * sres[i];
}
// magnetism
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = nrxx; i < 2 * nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta_mag * sres[i];
}
};
for (int i = nrxx; i < 2 * nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta_mag * sres[i];
}
};
this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true);
}
else if (GlobalV::NSPIN == 4)
{
// I do not merge this part with nspin=2 because the method of nspin=2 is almost done,
// while nspin=4 is not finished yet. I will try more methods for nspin=4 in the future.
rhor_in = chr->rho_save[0];
rhor_out = chr->rho[0];
const int nrxx = this->rhopw->nrxx;
auto screen = std::bind(&Charge_Mixing::Kerker_screen_real, this, std::placeholders::_1);
auto twobeta_mix
= [this, nrxx](double* out, const double* in, const double* sres) {
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = 0; i < nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta * sres[i];
}
// magnetism, mx, my, mz
#ifdef _OPENMP
#pragma omp parallel for schedule(static, 256)
#endif
for (int i = nrxx; i < 4 * nrxx; ++i)
{
out[i] = in[i] + this->mixing_beta_mag * sres[i];
}
};
this->mixing->push_data(this->rho_mdata, rhor_in, rhor_out, screen, twobeta_mix, true);
}

Expand Down Expand Up @@ -642,14 +701,18 @@ void Charge_Mixing::Kerker_screen_recip_new(std::complex<double>* drhog)
for (int is = 0; is < GlobalV::NSPIN; ++is)
{
// new mixing method only support nspin=2 not nspin=4
if (is == 1 && GlobalV::NSPIN == 2)
if (is >= 1)
{
if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
{
for (int ig = 0; ig < this->rhopw->npw; ig++)
{
drhog[is * this->rhopw->npw + ig] *= 1;
}
#ifdef __DEBUG
assert(is == 1); // make sure break works
#endif
double is_mag = GlobalV::NSPIN - 1;
//for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
//{
// drhog[is * this->rhopw->npw + ig] *= 1;
//}
break;
}
fac = GlobalV::MIXING_GG0_MAG;
Expand Down Expand Up @@ -692,13 +755,17 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
for (int is = 0; is < GlobalV::NSPIN; is++)
{

if (is == 1 && GlobalV::NSPIN == 2)
if (is >= 1)
{
if (GlobalV::MIXING_GG0_MAG <= 0.0001 || GlobalV::MIXING_BETA_MAG <= 0.1)
{
for (int ig = 0; ig < this->rhopw->npw; ig++)
#ifdef __DEBUG
assert(is == 1); // make sure break works
#endif
double is_mag = GlobalV::NSPIN - 1;
for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
{
drhog[is * this->rhopw->npw + ig] *= 0;
drhog[is * this->rhopw->npw + ig] = 0;
}
break;
}
Expand Down
22 changes: 13 additions & 9 deletions source/module_io/input.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -301,7 +301,7 @@ void Input::Default(void)
mixing_beta = -10;
mixing_ndim = 8;
mixing_gg0 = 1.00; // use Kerker defaultly
mixing_beta_mag = -10.0; // only set when nspin == 2
mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4
mixing_gg0_mag = 0.0; // defaultly exclude Kerker from mixing magnetic density
mixing_gg0_min = 0.1; // defaultly minimum kerker coefficient
mixing_tau = false;
Expand Down Expand Up @@ -2975,6 +2975,13 @@ void Input::Default_2(void) // jiyy add 2019-08-04
scf_thr_type = 1;
}
}

// set nspin with noncolin
if (noncolin || lspinorb)
{
nspin = 4;
}

// mixing parameters
if (mixing_beta < 0.0)
{
Expand All @@ -2990,28 +2997,25 @@ void Input::Default_2(void) // jiyy add 2019-08-04
}
else if (nspin == 4) // I will add this
{
mixing_beta = 0.2;
mixing_beta = 0.4;
mixing_beta_mag = 1.6;
mixing_gg0_mag = 0.0;
}
}
else
{
if (nspin == 2 && mixing_beta_mag < 0.0)
if ((nspin == 2 || nspin == 4) && mixing_beta_mag < 0.0)
{
if (mixing_beta <= 0.4)
{
mixing_beta_mag = 4 * mixing_beta;
}
else
{
mixing_beta_mag = 1.6;
mixing_beta_mag = 1.6; // 1.6 can be discussed
}
}
}
// set nspin with noncolin
if (noncolin || lspinorb)
{
nspin = 4;
}
}
#ifdef __MPI
void Input::Bcast()
Expand Down
6 changes: 3 additions & 3 deletions tests/integrate/104_PW_NC_magnetic/result.ref
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
etotref -5866.174730432665
etotperatomref -2933.0873652163
etotref -5866.197307043378
etotperatomref -2933.0986535217
pointgroupref O_h
spacegroupref O_h
nksibzref 1
totaltimeref 0.77367
totaltimeref 2.23
1 change: 1 addition & 0 deletions tests/integrate/204_NO_KP_NC_deltaspin/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ smearing_sigma 0.07

mixing_type pulay
mixing_beta 0.3
mixing_beta_mag 0.3
mixing_gg0 0.0

ks_solver genelpa
Expand Down