Skip to content

Commit 3b3466e

Browse files
sunliang98sunliang98
authored
Fix: Fix the Ewald force and stress when atom number of some elements are zero. (#5721)
* Fix: Fix the Ewald force when atom number is zero. * Fix: Fix the Ewald stress when atom number is zero. * Test: Add an integrate test 123_PW_zero_atom
1 parent 7b33110 commit 3b3466e

File tree

10 files changed

+186
-95
lines changed

10 files changed

+186
-95
lines changed

‎source/module_cell/atom_spec.cpp‎

Lines changed: 46 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -102,51 +102,54 @@ void Atom::bcast_atom(void)
102102
Parallel_Common::bcast_bool(this->flag_empty_element);
103103
Parallel_Common::bcast_double(mass);
104104

105-
if (GlobalV::MY_RANK != 0)
105+
if (na > 0)
106106
{
107-
assert(na != 0);
108-
this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
109-
this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
110-
this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
111-
this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
112-
this->mag.resize(na, 0);
113-
this->angle1.resize(na, 0);
114-
this->angle2.resize(na, 0);
115-
this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
116-
this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
117-
this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
118-
this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
119-
}
107+
if (GlobalV::MY_RANK != 0)
108+
{
109+
assert(na != 0);
110+
this->tau.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
111+
this->dis.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
112+
this->taud.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
113+
this->vel.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
114+
this->mag.resize(na, 0);
115+
this->angle1.resize(na, 0);
116+
this->angle2.resize(na, 0);
117+
this->m_loc_.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
118+
this->mbl.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
119+
this->lambda.resize(na, ModuleBase::Vector3<double>(0, 0, 0));
120+
this->constrain.resize(na, ModuleBase::Vector3<int>(0, 0, 0));
121+
}
120122

121-
for (int i = 0; i < na; i++)
122-
{
123-
Parallel_Common::bcast_double(tau[i].x);
124-
Parallel_Common::bcast_double(tau[i].y);
125-
Parallel_Common::bcast_double(tau[i].z);
126-
Parallel_Common::bcast_double(taud[i].x);
127-
Parallel_Common::bcast_double(taud[i].y);
128-
Parallel_Common::bcast_double(taud[i].z);
129-
Parallel_Common::bcast_double(dis[i].x);
130-
Parallel_Common::bcast_double(dis[i].y);
131-
Parallel_Common::bcast_double(dis[i].z);
132-
Parallel_Common::bcast_double(vel[i].x);
133-
Parallel_Common::bcast_double(vel[i].y);
134-
Parallel_Common::bcast_double(vel[i].z);
135-
Parallel_Common::bcast_double(mag[i]);
136-
Parallel_Common::bcast_double(angle1[i]);
137-
Parallel_Common::bcast_double(angle2[i]);
138-
Parallel_Common::bcast_double(m_loc_[i].x);
139-
Parallel_Common::bcast_double(m_loc_[i].y);
140-
Parallel_Common::bcast_double(m_loc_[i].z);
141-
Parallel_Common::bcast_int(mbl[i].x);
142-
Parallel_Common::bcast_int(mbl[i].y);
143-
Parallel_Common::bcast_int(mbl[i].z);
144-
Parallel_Common::bcast_double(lambda[i].x);
145-
Parallel_Common::bcast_double(lambda[i].y);
146-
Parallel_Common::bcast_double(lambda[i].z);
147-
Parallel_Common::bcast_int(constrain[i].x);
148-
Parallel_Common::bcast_int(constrain[i].y);
149-
Parallel_Common::bcast_int(constrain[i].z);
123+
for (int i = 0; i < na; i++)
124+
{
125+
Parallel_Common::bcast_double(tau[i].x);
126+
Parallel_Common::bcast_double(tau[i].y);
127+
Parallel_Common::bcast_double(tau[i].z);
128+
Parallel_Common::bcast_double(taud[i].x);
129+
Parallel_Common::bcast_double(taud[i].y);
130+
Parallel_Common::bcast_double(taud[i].z);
131+
Parallel_Common::bcast_double(dis[i].x);
132+
Parallel_Common::bcast_double(dis[i].y);
133+
Parallel_Common::bcast_double(dis[i].z);
134+
Parallel_Common::bcast_double(vel[i].x);
135+
Parallel_Common::bcast_double(vel[i].y);
136+
Parallel_Common::bcast_double(vel[i].z);
137+
Parallel_Common::bcast_double(mag[i]);
138+
Parallel_Common::bcast_double(angle1[i]);
139+
Parallel_Common::bcast_double(angle2[i]);
140+
Parallel_Common::bcast_double(m_loc_[i].x);
141+
Parallel_Common::bcast_double(m_loc_[i].y);
142+
Parallel_Common::bcast_double(m_loc_[i].z);
143+
Parallel_Common::bcast_int(mbl[i].x);
144+
Parallel_Common::bcast_int(mbl[i].y);
145+
Parallel_Common::bcast_int(mbl[i].z);
146+
Parallel_Common::bcast_double(lambda[i].x);
147+
Parallel_Common::bcast_double(lambda[i].y);
148+
Parallel_Common::bcast_double(lambda[i].z);
149+
Parallel_Common::bcast_int(constrain[i].x);
150+
Parallel_Common::bcast_int(constrain[i].y);
151+
Parallel_Common::bcast_int(constrain[i].z);
152+
}
150153
}
151154

152155
return;

‎source/module_cell/read_atoms.cpp‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -504,7 +504,18 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
504504
ModuleBase::WARNING("read_atom_positions", " atom number < 0.");
505505
return false;
506506
}
507-
if (na > 0)
507+
else if (na == 0)
508+
{
509+
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
510+
std::cout << " Warning: atom number is 0 for atom type: " << atoms[it].label << std::endl;
511+
std::cout << " If you are confident that this is not a mistake, please ignore this warning." << std::endl;
512+
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
513+
ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
514+
ofs_running << " Warning: atom number is 0 for atom type: " << atoms[it].label << std::endl;
515+
ofs_running << " If you are confident that this is not a mistake, please ignore this warning." << std::endl;
516+
ofs_running << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl;
517+
}
518+
else if (na > 0)
508519
{
509520
atoms[it].tau.resize(na, ModuleBase::Vector3<double>(0,0,0));
510521
atoms[it].dis.resize(na, ModuleBase::Vector3<double>(0,0,0));
@@ -891,6 +902,12 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
891902
ofs_running << std::endl;
892903
ModuleBase::GlobalFunc::OUT(ofs_running,"TOTAL ATOM NUMBER",nat);
893904

905+
if (nat == 0)
906+
{
907+
ModuleBase::WARNING("read_atom_positions","no atom in the system!");
908+
return false;
909+
}
910+
894911
// mohan add 2010-06-30
895912
//xiaohui modify 2015-03-15, cancel outputfile "STRU_READIN.xyz"
896913
//this->print_cell_xyz("STRU_READIN.xyz");

‎source/module_hamilt_pw/hamilt_pwdft/forces.cpp‎

Lines changed: 41 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -595,21 +595,24 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
595595
}
596596
for (int it = 0; it < ucell.ntype; it++)
597597
{
598-
double dzv;
599-
if (PARAM.inp.use_paw)
598+
if (ucell.atoms[it].na != 0)
600599
{
601-
#ifdef USE_PAW
602-
dzv = GlobalC::paw_cell.get_val(it);
603-
#endif
604-
}
605-
else
606-
{
607-
dzv = ucell.atoms[it].ncpp.zv;
608-
}
600+
double dzv;
601+
if (PARAM.inp.use_paw)
602+
{
603+
#ifdef USE_PAW
604+
dzv = GlobalC::paw_cell.get_val(it);
605+
#endif
606+
}
607+
else
608+
{
609+
dzv = ucell.atoms[it].ncpp.zv;
610+
}
609611

610-
for (int ig = igb; ig < ig_end; ++ig)
611-
{ // accumulate aux
612-
aux[ig] += dzv * conj(p_sf->strucFac(it, ig));
612+
for (int ig = igb; ig < ig_end; ++ig)
613+
{ // accumulate aux
614+
aux[ig] += dzv * conj(p_sf->strucFac(it, ig));
615+
}
613616
}
614617
}
615618
}
@@ -714,30 +717,33 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
714717
last_it = it;
715718
}
716719

717-
const auto ig_loop = [&](int ig_beg, int ig_end) {
718-
for (int ig = ig_beg; ig < ig_end; ig++)
719-
{
720-
const ModuleBase::Vector3<double> gcar = rho_basis->gcar[ig];
721-
const double arg = ModuleBase::TWO_PI * (gcar * ucell.atoms[it].tau[ia]);
722-
double sinp, cosp;
723-
ModuleBase::libm::sincos(arg, &sinp, &cosp);
724-
double sumnb = -cosp * aux[ig].imag() + sinp * aux[ig].real();
725-
forceion(iat, 0) += gcar[0] * sumnb;
726-
forceion(iat, 1) += gcar[1] * sumnb;
727-
forceion(iat, 2) += gcar[2] * sumnb;
728-
}
729-
};
720+
if (ucell.atoms[it].na != 0)
721+
{
722+
const auto ig_loop = [&](int ig_beg, int ig_end) {
723+
for (int ig = ig_beg; ig < ig_end; ig++)
724+
{
725+
const ModuleBase::Vector3<double> gcar = rho_basis->gcar[ig];
726+
const double arg = ModuleBase::TWO_PI * (gcar * ucell.atoms[it].tau[ia]);
727+
double sinp, cosp;
728+
ModuleBase::libm::sincos(arg, &sinp, &cosp);
729+
double sumnb = -cosp * aux[ig].imag() + sinp * aux[ig].real();
730+
forceion(iat, 0) += gcar[0] * sumnb;
731+
forceion(iat, 1) += gcar[1] * sumnb;
732+
forceion(iat, 2) += gcar[2] * sumnb;
733+
}
734+
};
730735

731-
// skip ig_gge0 point by separating ig loop into two part
732-
ig_loop(0, ig_gap);
733-
ig_loop(ig_gap + 1, rho_basis->npw);
736+
// skip ig_gge0 point by separating ig loop into two part
737+
ig_loop(0, ig_gap);
738+
ig_loop(ig_gap + 1, rho_basis->npw);
734739

735-
forceion(iat, 0) *= it_fact;
736-
forceion(iat, 1) *= it_fact;
737-
forceion(iat, 2) *= it_fact;
740+
forceion(iat, 0) *= it_fact;
741+
forceion(iat, 1) *= it_fact;
742+
forceion(iat, 2) *= it_fact;
738743

739-
++iat;
740-
ucell.step_iait(&ia, &it);
744+
++iat;
745+
ucell.step_iait(&ia, &it);
746+
}
741747
}
742748

743749
// means that the processor contains G=0 term.
@@ -771,7 +777,7 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
771777
int T2 = 0;
772778
while (iat2 < this->nat)
773779
{
774-
if (iat1 != iat2)
780+
if (iat1 != iat2 && ucell.atoms[T2].na != 0 && ucell.atoms[T1].na != 0)
775781
{
776782
ModuleBase::Vector3<double> d_tau
777783
= ucell.atoms[T1].tau[I1] - ucell.atoms[T2].tau[I2];

‎source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp‎

Lines changed: 19 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -137,25 +137,28 @@ void Stress_Func<FPTYPE, Device>::stress_ewa(const UnitCell& ucell,
137137

138138
while (ijat < ijat_end)
139139
{
140-
//calculate tau[na]-tau[nb]
141-
d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j];
142-
//generates nearest-neighbors shells
143-
H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm);
144-
for(int nr=0 ; nr<nrm ; nr++)
140+
if (ucell.atoms[it].na != 0 && ucell.atoms[jt].na != 0)
145141
{
146-
rr=sqrt(r2[nr]) * ucell.lat0;
147-
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)));
148-
for(int l=0; l<3; l++)
142+
//calculate tau[na]-tau[nb]
143+
d_tau = ucell.atoms[it].tau[i] - ucell.atoms[jt].tau[j];
144+
//generates nearest-neighbors shells
145+
H_Ewald_pw::rgen(d_tau, rmax, irr, ucell.latvec, ucell.G, r, r2, nrm);
146+
for(int nr=0 ; nr<nrm ; nr++)
149147
{
150-
for(int m=0; m<l+1; m++)
148+
rr=sqrt(r2[nr]) * ucell.lat0;
149+
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)));
150+
for(int l=0; l<3; l++)
151151
{
152-
r0[0] = r[nr].x;
153-
r0[1] = r[nr].y;
154-
r0[2] = r[nr].z;
155-
local_sigma(l,m) += fac * r0[l] * r0[m];
156-
}//end m
157-
}//end l
158-
}//end nr
152+
for(int m=0; m<l+1; m++)
153+
{
154+
r0[0] = r[nr].x;
155+
r0[1] = r[nr].y;
156+
r0[2] = r[nr].z;
157+
local_sigma(l,m) += fac * r0[l] * r0[m];
158+
}//end m
159+
}//end l
160+
}//end nr
161+
}
159162

160163
++ijat;
161164
ucell.step_jajtiait(&j, &jt, &i, &it);

‎tests/integrate/123_PW_zero_atom/INPUT‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
INPUT_PARAMETERS
2+
#Parameters (1.General)
3+
suffix autotest
4+
calculation scf
5+
6+
nbands 6
7+
symmetry 0
8+
pseudo_dir ../../PP_ORB/
9+
10+
#Parameters (2.Iteration)
11+
ecutwfc 20
12+
13+
#Parameters (3.Basis)
14+
basis_type pw
15+
16+
#Parameters (4.Smearing)
17+
smearing_method gauss
18+
smearing_sigma 0.0002
19+
20+
#Parameters (5.Mixing)
21+
mixing_type broyden
22+
mixing_beta 0.7
23+
cal_force 1
24+
test_force 1
25+
cal_stress 1
26+
test_stress 1
27+
28+
pw_seed 1

‎tests/integrate/123_PW_zero_atom/KPT‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
K_POINTS
2+
0
3+
Gamma
4+
2 2 2 0 0 0

‎tests/integrate/123_PW_zero_atom/STRU‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
ATOMIC_SPECIES
2+
Al 13 Al_ONCV_PBE-1.0.upf upf201
3+
Si 14 Si_ONCV_PBE-1.0.upf upf201
4+
5+
LATTICE_CONSTANT
6+
10.2 // add lattice constant
7+
8+
LATTICE_VECTORS
9+
0.5 0.5 0.0
10+
0.5 0.0 0.5
11+
0.0 0.5 0.5
12+
13+
ATOMIC_POSITIONS
14+
Direct
15+
16+
Al
17+
0.0
18+
0
19+
Si // Element type
20+
0.0 // magnetism
21+
2
22+
0.00 0.00 0.00 1 1 1
23+
0.25 0.25 0.25 1 1 1

‎tests/integrate/123_PW_zero_atom/jd‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
test calculation when atom number is zero.

‎tests/integrate/123_PW_zero_atom/result.ref‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
etotref -211.8003327929410489
2+
etotperatomref -105.9001663965
3+
totalforceref 0.000014
4+
totalstressref 368.726447
5+
totaltimeref 0.61

‎tests/integrate/CASES_CPU.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,7 @@
9898
121_PW_KPAR
9999
121_PW_kspacing
100100
127_PW_15_PK_AF
101+
123_PW_zero_atom
101102
128_PW_zero_ntype
102103
133_PW_DJ_PK
103104
135_PW_15_PK

0 commit comments

Comments
 (0)