Skip to content

Fix: Fix the Ewald force and stress when atom number of some elements are zero. #5721

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
mohanchen merged 4 commits into from
Dec 12, 2024
Merged

Fix: Fix the Ewald force and stress when atom number of some elements are zero. #5721

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
95c2b49
Fix: Fix the Ewald force when atom number is zero.
sunliang98 Dec 11, 2024
594833c
Fix: Fix the Ewald stress when atom number is zero.
sunliang98 Dec 11, 2024
be0f278
Test: Add an integrate test 123_PW_zero_atom
sunliang98 Dec 11, 2024
7dbf6f4
Merge branch 'develop' of https://github.com/sunliang98/abacus-develo…
sunliang98 Dec 11, 2024
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
89 changes: 46 additions & 43 deletions source/module_cell/atom_spec.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,51 +102,54 @@ void Atom::bcast_atom(void)
Parallel_Common::bcast_bool(this->flag_empty_element);
Parallel_Common::bcast_double(mass);

if (GlobalV::MY_RANK != 0)
if (na > 0)
{
assert(na != 0);
this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->mag.resize(na, 0);
this->angle1.resize(na, 0);
this->angle2.resize(na, 0);
this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
}
if (GlobalV::MY_RANK != 0)
{
assert(na != 0);
this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->mag.resize(na, 0);
this->angle1.resize(na, 0);
this->angle2.resize(na, 0);
this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
}

for (int i = 0; i < na; i++)
{
Parallel_Common::bcast_double(tau[i].x);
Parallel_Common::bcast_double(tau[i].y);
Parallel_Common::bcast_double(tau[i].z);
Parallel_Common::bcast_double(taud[i].x);
Parallel_Common::bcast_double(taud[i].y);
Parallel_Common::bcast_double(taud[i].z);
Parallel_Common::bcast_double(dis[i].x);
Parallel_Common::bcast_double(dis[i].y);
Parallel_Common::bcast_double(dis[i].z);
Parallel_Common::bcast_double(vel[i].x);
Parallel_Common::bcast_double(vel[i].y);
Parallel_Common::bcast_double(vel[i].z);
Parallel_Common::bcast_double(mag[i]);
Parallel_Common::bcast_double(angle1[i]);
Parallel_Common::bcast_double(angle2[i]);
Parallel_Common::bcast_double(m_loc_[i].x);
Parallel_Common::bcast_double(m_loc_[i].y);
Parallel_Common::bcast_double(m_loc_[i].z);
Parallel_Common::bcast_int(mbl[i].x);
Parallel_Common::bcast_int(mbl[i].y);
Parallel_Common::bcast_int(mbl[i].z);
Parallel_Common::bcast_double(lambda[i].x);
Parallel_Common::bcast_double(lambda[i].y);
Parallel_Common::bcast_double(lambda[i].z);
Parallel_Common::bcast_int(constrain[i].x);
Parallel_Common::bcast_int(constrain[i].y);
Parallel_Common::bcast_int(constrain[i].z);
for (int i = 0; i < na; i++)
{
Parallel_Common::bcast_double(tau[i].x);
Parallel_Common::bcast_double(tau[i].y);
Parallel_Common::bcast_double(tau[i].z);
Parallel_Common::bcast_double(taud[i].x);
Parallel_Common::bcast_double(taud[i].y);
Parallel_Common::bcast_double(taud[i].z);
Parallel_Common::bcast_double(dis[i].x);
Parallel_Common::bcast_double(dis[i].y);
Parallel_Common::bcast_double(dis[i].z);
Parallel_Common::bcast_double(vel[i].x);
Parallel_Common::bcast_double(vel[i].y);
Parallel_Common::bcast_double(vel[i].z);
Parallel_Common::bcast_double(mag[i]);
Parallel_Common::bcast_double(angle1[i]);
Parallel_Common::bcast_double(angle2[i]);
Parallel_Common::bcast_double(m_loc_[i].x);
Parallel_Common::bcast_double(m_loc_[i].y);
Parallel_Common::bcast_double(m_loc_[i].z);
Parallel_Common::bcast_int(mbl[i].x);
Parallel_Common::bcast_int(mbl[i].y);
Parallel_Common::bcast_int(mbl[i].z);
Parallel_Common::bcast_double(lambda[i].x);
Parallel_Common::bcast_double(lambda[i].y);
Parallel_Common::bcast_double(lambda[i].z);
Parallel_Common::bcast_int(constrain[i].x);
Parallel_Common::bcast_int(constrain[i].y);
Parallel_Common::bcast_int(constrain[i].z);
}
}

return;
Expand Down
19 changes: 18 additions & 1 deletion source/module_cell/read_atoms.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -504,7 +504,18 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
ModuleBase::WARNING("read_atom_positions", " atom number < 0.");
return false;
}
if (na > 0)
else if (na == 0)
{
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
std::cout << " Warning: atom number is 0 for atom type: " << atoms[it].label << std::endl;
std::cout << " If you are confident that this is not a mistake, please ignore this warning." << std::endl;
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
ofs_running << " Warning: atom number is 0 for atom type: " << atoms[it].label << std::endl;
ofs_running << " If you are confident that this is not a mistake, please ignore this warning." << std::endl;
ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
}
else if (na > 0)
{
atoms[it].tau.resize(na, ModuleBase::Vector3<double>(0,0,0));
atoms[it].dis.resize(na, ModuleBase::Vector3<double>(0,0,0));
Expand Down Expand Up @@ -891,6 +902,12 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
ofs_running << std::endl;
ModuleBase::GlobalFunc::OUT(ofs_running,"TOTAL ATOM NUMBER",nat);

if (nat == 0)
{
ModuleBase::WARNING("read_atom_positions","no atom in the system!");
return false;
}

// mohan add 2010-06-30
//xiaohui modify 2015-03-15, cancel outputfile "STRU_READIN.xyz"
//this->print_cell_xyz("STRU_READIN.xyz");
Expand Down
76 changes: 41 additions & 35 deletions source/module_hamilt_pw/hamilt_pwdft/forces.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -595,21 +595,24 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
}
for (int it = 0; it < ucell.ntype; it++)
{
double dzv;
if (PARAM.inp.use_paw)
if (ucell.atoms[it].na != 0)
{
#ifdef USE_PAW
dzv = GlobalC::paw_cell.get_val(it);
#endif
}
else
{
dzv = ucell.atoms[it].ncpp.zv;
}
double dzv;
if (PARAM.inp.use_paw)
{
#ifdef USE_PAW
dzv = GlobalC::paw_cell.get_val(it);
#endif
}
else
{
dzv = ucell.atoms[it].ncpp.zv;
}

for (int ig = igb; ig < ig_end; ++ig)
{ // accumulate aux
aux[ig] += dzv * conj(p_sf->strucFac(it, ig));
for (int ig = igb; ig < ig_end; ++ig)
{ // accumulate aux
aux[ig] += dzv * conj(p_sf->strucFac(it, ig));
}
}
}
}
Expand Down Expand Up @@ -714,30 +717,33 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
last_it = it;
}

const auto ig_loop = [&](int ig_beg, int ig_end) {
for (int ig = ig_beg; ig < ig_end; ig++)
{
const ModuleBase::Vector3<double> gcar = rho_basis->gcar[ig];
const double arg = ModuleBase::TWO_PI * (gcar * ucell.atoms[it].tau[ia]);
double sinp, cosp;
ModuleBase::libm::sincos(arg, &sinp, &cosp);
double sumnb = -cosp * aux[ig].imag() + sinp * aux[ig].real();
forceion(iat, 0) += gcar[0] * sumnb;
forceion(iat, 1) += gcar[1] * sumnb;
forceion(iat, 2) += gcar[2] * sumnb;
}
};
if (ucell.atoms[it].na != 0)
{
const auto ig_loop = [&](int ig_beg, int ig_end) {
for (int ig = ig_beg; ig < ig_end; ig++)
{
const ModuleBase::Vector3<double> gcar = rho_basis->gcar[ig];
const double arg = ModuleBase::TWO_PI * (gcar * ucell.atoms[it].tau[ia]);
double sinp, cosp;
ModuleBase::libm::sincos(arg, &sinp, &cosp);
double sumnb = -cosp * aux[ig].imag() + sinp * aux[ig].real();
forceion(iat, 0) += gcar[0] * sumnb;
forceion(iat, 1) += gcar[1] * sumnb;
forceion(iat, 2) += gcar[2] * sumnb;
}
};

// skip ig_gge0 point by separating ig loop into two part
ig_loop(0, ig_gap);
ig_loop(ig_gap + 1, rho_basis->npw);
// skip ig_gge0 point by separating ig loop into two part
ig_loop(0, ig_gap);
ig_loop(ig_gap + 1, rho_basis->npw);

forceion(iat, 0) *= it_fact;
forceion(iat, 1) *= it_fact;
forceion(iat, 2) *= it_fact;
forceion(iat, 0) *= it_fact;
forceion(iat, 1) *= it_fact;
forceion(iat, 2) *= it_fact;

++iat;
ucell.step_iait(&ia, &it);
++iat;
ucell.step_iait(&ia, &it);
}
}

// means that the processor contains G=0 term.
Expand Down Expand Up @@ -771,7 +777,7 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
int T2 = 0;
while (iat2 < this->nat)
{
if (iat1 != iat2)
if (iat1 != iat2 && ucell.atoms[T2].na != 0 && ucell.atoms[T1].na != 0)
{
ModuleBase::Vector3<double> d_tau
= ucell.atoms[T1].tau[I1] - ucell.atoms[T2].tau[I2];
Expand Down
35 changes: 19 additions & 16 deletions source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,25 +137,28 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(const UnitCell& ucell,

while (ijat < ijat_end)
{
//calculate tau[na]-tau[nb]
d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j];
//generates nearest-neighbors shells
H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm);
for(int nr=0 ; nr<nrm ; nr++)
if (ucell.atoms[it].na != 0 && ucell.atoms[jt].na != 0)
{
rr=sqrt(r2[nr]) * ucell.lat0;
fac = -ModuleBase::e2/2.0/ucell.omega*pow(ucell.lat0,2)*ucell.atoms[it].ncpp.zv * ucell.atoms[jt].ncpp.zv / pow(rr,3) * (erfc(sqa*rr)+rr * sq8a_2pi * ModuleBase::libm::exp(-alpha * pow(rr,2)));
for(int l=0; l<3; l++)
//calculate tau[na]-tau[nb]
d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j];
//generates nearest-neighbors shells
H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm);
for(int nr=0 ; nr<nrm ; nr++)
{
for(int m=0; m<l+1; m++)
rr=sqrt(r2[nr]) * ucell.lat0;
fac = -ModuleBase::e2/2.0/ucell.omega*pow(ucell.lat0,2)*ucell.atoms[it].ncpp.zv * ucell.atoms[jt].ncpp.zv / pow(rr,3) * (erfc(sqa*rr)+rr * sq8a_2pi * ModuleBase::libm::exp(-alpha * pow(rr,2)));
for(int l=0; l<3; l++)
{
r0[0] = r[nr].x;
r0[1] = r[nr].y;
r0[2] = r[nr].z;
local_sigma(l,m) += fac * r0[l] * r0[m];
}//end m
}//end l
}//end nr
for(int m=0; m<l+1; m++)
{
r0[0] = r[nr].x;
r0[1] = r[nr].y;
r0[2] = r[nr].z;
local_sigma(l,m) += fac * r0[l] * r0[m];
}//end m
}//end l
}//end nr
}

++ijat;
ucell.step_jajtiait(&j, &jt, &i, &it);
Expand Down
28 changes: 28 additions & 0 deletions tests/integrate/123_PW_zero_atom/INPUT
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
INPUT_PARAMETERS
#Parameters (1.General)
suffix autotest
calculation scf

nbands 6
symmetry 0
pseudo_dir ../../PP_ORB/

#Parameters (2.Iteration)
ecutwfc 20

#Parameters (3.Basis)
basis_type pw

#Parameters (4.Smearing)
smearing_method gauss
smearing_sigma 0.0002

#Parameters (5.Mixing)
mixing_type broyden
mixing_beta 0.7
cal_force 1
test_force 1
cal_stress 1
test_stress 1

pw_seed 1
4 changes: 4 additions & 0 deletions tests/integrate/123_PW_zero_atom/KPT
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
K_POINTS
0
Gamma
2 2 2 0 0 0
23 changes: 23 additions & 0 deletions tests/integrate/123_PW_zero_atom/STRU
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
ATOMIC_SPECIES
Al 13 Al_ONCV_PBE-1.0.upf upf201
Si 14 Si_ONCV_PBE-1.0.upf upf201

LATTICE_CONSTANT
10.2 // add lattice constant

LATTICE_VECTORS
0.5 0.5 0.0
0.5 0.0 0.5
0.0 0.5 0.5

ATOMIC_POSITIONS
Direct

Al
0.0
0
Si // Element type
0.0 // magnetism
2
0.00 0.00 0.00 1 1 1
0.25 0.25 0.25 1 1 1
1 change: 1 addition & 0 deletions tests/integrate/123_PW_zero_atom/jd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
test calculation when atom number is zero.
5 changes: 5 additions & 0 deletions tests/integrate/123_PW_zero_atom/result.ref
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
etotref -211.8003327929410489
etotperatomref -105.9001663965
totalforceref 0.000014
totalstressref 368.726447
totaltimeref 0.61
1 change: 1 addition & 0 deletions tests/integrate/CASES_CPU.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@
121_PW_KPAR
121_PW_kspacing
127_PW_15_PK_AF
123_PW_zero_atom
128_PW_zero_ntype
133_PW_DJ_PK
135_PW_15_PK
Expand Down