From 95c2b49f528f3b7cbfdb89d81be89a35f184e4b0 Mon Sep 17 00:00:00 2001 From: sunliang98 <1700011430@pku.edu.cn> Date: Wed, 11 Dec 2024 20:09:58 +0800 Subject: [PATCH 1/3] Fix: Fix the Ewald force when atom number is zero. --- source/module_cell/atom_spec.cpp | 89 ++++++++++--------- source/module_cell/read_atoms.cpp | 19 +++- .../module_hamilt_pw/hamilt_pwdft/forces.cpp | 76 ++++++++-------- 3 files changed, 105 insertions(+), 79 deletions(-) diff --git a/source/module_cell/atom_spec.cpp b/source/module_cell/atom_spec.cpp index f00a062c1d..0bf5045f2d 100644 --- a/source/module_cell/atom_spec.cpp +++ b/source/module_cell/atom_spec.cpp @@ -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(0, 0, 0)); - this->dis.resize(na, ModuleBase::Vector3(0, 0, 0)); - this->taud.resize(na, ModuleBase::Vector3(0, 0, 0)); - this->vel.resize(na, ModuleBase::Vector3(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(0, 0, 0)); - this->mbl.resize(na, ModuleBase::Vector3(0, 0, 0)); - this->lambda.resize(na, ModuleBase::Vector3(0, 0, 0)); - this->constrain.resize(na, ModuleBase::Vector3(0, 0, 0)); - } + if (GlobalV::MY_RANK != 0) + { + assert(na != 0); + this->tau.resize(na, ModuleBase::Vector3(0, 0, 0)); + this->dis.resize(na, ModuleBase::Vector3(0, 0, 0)); + this->taud.resize(na, ModuleBase::Vector3(0, 0, 0)); + this->vel.resize(na, ModuleBase::Vector3(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(0, 0, 0)); + this->mbl.resize(na, ModuleBase::Vector3(0, 0, 0)); + this->lambda.resize(na, ModuleBase::Vector3(0, 0, 0)); + this->constrain.resize(na, ModuleBase::Vector3(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; diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index 81608ce609..30ae5c504b 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -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(0,0,0)); atoms[it].dis.resize(na, ModuleBase::Vector3(0,0,0)); @@ -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"); diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp index c1fcd2299c..dd24ce82c7 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp @@ -595,21 +595,24 @@ void Forces::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)); + } } } } @@ -714,30 +717,33 @@ void Forces::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 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 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. @@ -771,7 +777,7 @@ void Forces::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 d_tau = ucell.atoms[T1].tau[I1] - ucell.atoms[T2].tau[I2]; From 594833cbb4064e16aa6affc2ac7353f143278467 Mon Sep 17 00:00:00 2001 From: sunliang98 <1700011430@pku.edu.cn> Date: Wed, 11 Dec 2024 20:20:21 +0800 Subject: [PATCH 2/3] Fix: Fix the Ewald stress when atom number is zero. --- .../hamilt_pwdft/stress_func_ewa.cpp | 35 ++++++++++--------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp index fd6bf8dc0a..5dbe2cd6b4 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_ewa.cpp @@ -137,25 +137,28 @@ void Stress_Func::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 Date: Wed, 11 Dec 2024 20:45:05 +0800 Subject: [PATCH 3/3] Test: Add an integrate test 123_PW_zero_atom --- tests/integrate/123_PW_zero_atom/INPUT | 28 +++++++++++++++++++++ tests/integrate/123_PW_zero_atom/KPT | 4 +++ tests/integrate/123_PW_zero_atom/STRU | 23 +++++++++++++++++ tests/integrate/123_PW_zero_atom/jd | 1 + tests/integrate/123_PW_zero_atom/result.ref | 5 ++++ tests/integrate/CASES_CPU.txt | 1 + 6 files changed, 62 insertions(+) create mode 100644 tests/integrate/123_PW_zero_atom/INPUT create mode 100644 tests/integrate/123_PW_zero_atom/KPT create mode 100644 tests/integrate/123_PW_zero_atom/STRU create mode 100644 tests/integrate/123_PW_zero_atom/jd create mode 100644 tests/integrate/123_PW_zero_atom/result.ref diff --git a/tests/integrate/123_PW_zero_atom/INPUT b/tests/integrate/123_PW_zero_atom/INPUT new file mode 100644 index 0000000000..8488f61585 --- /dev/null +++ b/tests/integrate/123_PW_zero_atom/INPUT @@ -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 \ No newline at end of file diff --git a/tests/integrate/123_PW_zero_atom/KPT b/tests/integrate/123_PW_zero_atom/KPT new file mode 100644 index 0000000000..f5f7f4ec34 --- /dev/null +++ b/tests/integrate/123_PW_zero_atom/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 2 0 0 0 diff --git a/tests/integrate/123_PW_zero_atom/STRU b/tests/integrate/123_PW_zero_atom/STRU new file mode 100644 index 0000000000..49a2ef9cb7 --- /dev/null +++ b/tests/integrate/123_PW_zero_atom/STRU @@ -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 diff --git a/tests/integrate/123_PW_zero_atom/jd b/tests/integrate/123_PW_zero_atom/jd new file mode 100644 index 0000000000..ebecbd81ec --- /dev/null +++ b/tests/integrate/123_PW_zero_atom/jd @@ -0,0 +1 @@ +test calculation when atom number is zero. diff --git a/tests/integrate/123_PW_zero_atom/result.ref b/tests/integrate/123_PW_zero_atom/result.ref new file mode 100644 index 0000000000..c93cf82f0f --- /dev/null +++ b/tests/integrate/123_PW_zero_atom/result.ref @@ -0,0 +1,5 @@ +etotref -211.8003327929410489 +etotperatomref -105.9001663965 +totalforceref 0.000014 +totalstressref 368.726447 +totaltimeref 0.61 diff --git a/tests/integrate/CASES_CPU.txt b/tests/integrate/CASES_CPU.txt index 8277d51eab..1c76104ae0 100644 --- a/tests/integrate/CASES_CPU.txt +++ b/tests/integrate/CASES_CPU.txt @@ -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