From 0d17a87095f02a7e9069e3e37cec2c32f17670b7 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 3 Mar 2025 12:15:24 +0800 Subject: [PATCH 01/14] first version --- source/module_cell/check_atomic_stru.cpp | 122 +++++++++++++---------- 1 file changed, 71 insertions(+), 51 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 527cd00e77..280beca780 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -1,5 +1,5 @@ #include "check_atomic_stru.h" - +#include "module_base/timer.h" #include "module_base/element_covalent_radius.h" void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { @@ -7,65 +7,73 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) // and compare with the covalent_bond_length, // if there has bond length is shorter than covalent_bond_length * factor, // we think this structure is unreasonable. - const double warning_coef = 0.6; + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); assert(ucell.ntype > 0); - std::stringstream errorlog; + const int ntype = ucell.ntype; + const double lat0 = ucell.lat0; bool all_pass = true; bool no_warning = true; - for (int it1 = 0; it1 < ucell.ntype; it1++) { + const double warning_coef = 0.6; + const double min_factor_coef=std::min(warning_coef,factor); + + std::stringstream errorlog; + errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); + + std::vector symbol_covalent_radiuss(ntype); + for (int it = 0;it < ntype ; it++) + { std::string symbol1 = ""; - for (char ch: ucell.atoms[it1].label) { - if (std::isalpha(ch)) { + for (char ch: ucell.atoms[it].label) + { + if (std::isalpha(ch)) + { symbol1.push_back(ch); } } - // std::string symbol1 = ucell.atoms[it1].label; - double symbol1_covalent_radius; - if (ModuleBase::CovalentRadius.find(symbol1) - != ModuleBase::CovalentRadius.end()) { - symbol1_covalent_radius = ModuleBase::CovalentRadius.at(symbol1); - } else { + + if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) + { + symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); + } + else + { std::stringstream mess; - mess << "Notice: symbol '" << symbol1 - << "' is not an element symbol!!!! "; + mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; mess << "set the covalent radius to be 0." << std::endl; GlobalV::ofs_running << mess.str(); std::cout << mess.str(); - symbol1_covalent_radius = 0.0; } + } + for (int it1 = 0; it1 < ntype; it1++) { + const double symbol1_covalent_radius=symbol_covalent_radiuss[it1]; for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) { double x1 = ucell.atoms[it1].taud[ia1].x; double y1 = ucell.atoms[it1].taud[ia1].y; double z1 = ucell.atoms[it1].taud[ia1].z; - for (int it2 = 0; it2 < ucell.ntype; it2++) { + for (int it2 = it1; it2 < ucell.ntype; it2++) { std::string symbol2 = ucell.atoms[it2].label; - double symbol2_covalent_radius; - if (ModuleBase::CovalentRadius.find(symbol2) - != ModuleBase::CovalentRadius.end()) { - symbol2_covalent_radius - = ModuleBase::CovalentRadius.at(symbol2); - } else { - symbol2_covalent_radius = 0.0; - } - - double covalent_length - = (symbol1_covalent_radius + symbol2_covalent_radius) - / ModuleBase::BOHR_TO_A; - - for (int ia2 = 0; ia2 < ucell.atoms[it2].na; ia2++) { + double symbol2_covalent_radius=symbol_covalent_radiuss[it1]; + + double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / ModuleBase::BOHR_TO_A; + const double min_error = covalent_length * min_factor_coef; + + + for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) + { + double delta_x = ucell.atoms[it2].taud[ia2].x - x1; + double delta_y = ucell.atoms[it2].taud[ia2].y - y1; + double delta_z = ucell.atoms[it2].taud[ia2].z - z1; + for (int a = -1; a < 2; a++) { for (int b = -1; b < 2; b++) { for (int c = -1; c < 2; c++) { - if (it1 > it2) { - continue; - } else if (it1 == it2 && ia1 > ia2) { - continue; - } else if (it1 == it2 && ia1 == ia2 && a == 0 - && b == 0 && c == 0) { + if (it1 == it2 && ia1 == ia2 && a == 0 + && b == 0 && c == 0) + { continue; -} + } double x2 = ucell.atoms[it2].taud[ia2].x + a; double y2 = ucell.atoms[it2].taud[ia2].y + b; @@ -86,11 +94,8 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) 2)) * ucell.lat0; - if (bond_length < covalent_length * factor - || bond_length - < covalent_length * warning_coef) { - errorlog.setf(std::ios_base::fixed, - std::ios_base::floatfield); + if (bond_length < min_error) { + errorlog << std::setw(3) << ia1 + 1 << "-th " << std::setw(3) << ucell.atoms[it1].label << ", "; @@ -105,14 +110,17 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) << bond_length << " Bohr ("; errorlog << bond_length * ModuleBase::BOHR_TO_A - << " Angstrom)" << std::endl; - - if (bond_length - < covalent_length * factor) { - all_pass = false; - } else { - no_warning = false; + << " Angstrom)\n"; + if (all_pass || no_warning) + { + if (bond_length < covalent_length * factor) + { + all_pass = false; + } else { + no_warning = false; + } } + } } // c } // b @@ -121,7 +129,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) } // it2 } // ia1 } // it1 - + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); if (!all_pass || !no_warning) { std::stringstream mess; mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" @@ -143,7 +151,18 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); std::cout << mess.str() << mess.str() << mess.str() << std::endl; + std::ofstream outFile("example2.txt"); + if (outFile.is_open()) { + + // 将stringstream中的内容写入到文件 + outFile << errorlog.str(); + // 关闭文件 + outFile.close(); + std::cout << "数据成功写入到文件。\n"; + } else { + std::cerr << "无法打开文件进行写入。\n"; + } if (!all_pass) { mess.clear(); mess.str(""); @@ -154,7 +173,8 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) << ") in INPUT file." << std::endl; GlobalV::ofs_running << mess.str(); std::cout << mess.str(); - ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); + // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } } + } From cfd37789c9f9e174cedefecb4f31694e8f4988f0 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 3 Mar 2025 13:24:22 +0800 Subject: [PATCH 02/14] change add --- source/module_cell/check_atomic_stru.cpp | 116 +++++++++++------------ 1 file changed, 56 insertions(+), 60 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 280beca780..d4c679db23 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -43,8 +43,21 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) GlobalV::ofs_running << mess.str(); std::cout << mess.str(); } + } + double * latvec = new double [9]; + latvec[0]=ucell.a1.x; latvec[1]=ucell.a2.x; latvec[2]=ucell.a3.x; + latvec[3]=ucell.a1.y; latvec[4]=ucell.a2.y; latvec[5]=ucell.a3.y; + latvec[6]=ucell.a1.z; latvec[7]=ucell.a2.z; latvec[8]=ucell.a3.z; + std::vector A(27*3); + for (int i=0;i<27;i++) + { + int a = (i / 9) % 3 - 1; // 计算a的值 + int b = (i / 3) % 3 - 1; // 计算b的值 + int c = i % 3 - 1; // 计算c的值 + A[3*i]= a*latvec[0] + b *latvec[1] + c* latvec[2]; + A[3*i+1]= a*latvec[3] + b *latvec[4] + c* latvec[5]; + A[3*i+2]= a*latvec[6] + b *latvec[7] + c* latvec[8]; } - for (int it1 = 0; it1 < ntype; it1++) { const double symbol1_covalent_radius=symbol_covalent_radiuss[it1]; for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) { @@ -57,78 +70,61 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) double symbol2_covalent_radius=symbol_covalent_radiuss[it1]; double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / ModuleBase::BOHR_TO_A; - const double min_error = covalent_length * min_factor_coef; - + const double min_error = covalent_length * min_factor_coef / ucell.lat0; + const double min_error_2 = min_error * min_error; for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) { double delta_x = ucell.atoms[it2].taud[ia2].x - x1; double delta_y = ucell.atoms[it2].taud[ia2].y - y1; double delta_z = ucell.atoms[it2].taud[ia2].z - z1; - - for (int a = -1; a < 2; a++) { - for (int b = -1; b < 2; b++) { - for (int c = -1; c < 2; c++) { - if (it1 == it2 && ia1 == ia2 && a == 0 - && b == 0 && c == 0) + const double delta_x_lat = delta_x*latvec[0] + delta_y*latvec[1]+ delta_z*latvec[2]; + const double delta_y_lat = delta_x*latvec[3] + delta_y*latvec[4]+ delta_z*latvec[5]; + const double delta_z_lat = delta_x*latvec[6] + delta_y*latvec[7]+ delta_z*latvec[8]; + for (int i=0;i<27;i++) + { + + if (it1 == it2 && ia1 == ia2 && i==13) + { + continue; + } + double part1= delta_x_lat+A[3*i]; + double part2= delta_y_lat+A[3*i+1]; + double part3= delta_z_lat+A[3*i+2]; + double bond_length = part1*part1 + part2*part2+ part3*part3; + + if (bond_length < min_error_2) + { + errorlog << std::setw(3) << ia1 + 1 + << "-th " << std::setw(3) + << ucell.atoms[it1].label << ", " + << std::setw(3) << ia2 + 1 + << "-th " << std::setw(3) + << ucell.atoms[it2].label + << " (cell:" << std::setw(2) << (i / 9) % 3 - 1 + << " " << std::setw(2) << (i / 3) % 3 - 1 << " " + << std::setw(2) << i % 3 - 1 << "), distance= " + << std::setprecision(3) + << sqrt(bond_length)/lat0 << " Bohr (" + << sqrt(bond_length) /lat0 * ModuleBase::BOHR_TO_A + << " Angstrom)\n"; + if (all_pass || no_warning) + { + if (bond_length < pow(covalent_length * factor/lat0,2)) { - continue; - } - - double x2 = ucell.atoms[it2].taud[ia2].x + a; - double y2 = ucell.atoms[it2].taud[ia2].y + b; - double z2 = ucell.atoms[it2].taud[ia2].z + c; - - double bond_length - = sqrt(pow((x2 - x1) * ucell.a1.x - + (y2 - y1) * ucell.a2.x - + (z2 - z1) * ucell.a3.x, - 2) - + pow((x2 - x1) * ucell.a1.y - + (y2 - y1) * ucell.a2.y - + (z2 - z1) * ucell.a3.y, - 2) - + pow((x2 - x1) * ucell.a1.z - + (y2 - y1) * ucell.a2.z - + (z2 - z1) * ucell.a3.z, - 2)) - * ucell.lat0; - - if (bond_length < min_error) { - - errorlog << std::setw(3) << ia1 + 1 - << "-th " << std::setw(3) - << ucell.atoms[it1].label << ", "; - errorlog << std::setw(3) << ia2 + 1 - << "-th " << std::setw(3) - << ucell.atoms[it2].label; - errorlog << " (cell:" << std::setw(2) << a - << " " << std::setw(2) << b << " " - << std::setw(2) << c << ")"; - errorlog << ", distance= " - << std::setprecision(3) - << bond_length << " Bohr ("; - errorlog - << bond_length * ModuleBase::BOHR_TO_A - << " Angstrom)\n"; - if (all_pass || no_warning) - { - if (bond_length < covalent_length * factor) - { - all_pass = false; - } else { - no_warning = false; - } - } - + all_pass = false; + } else { + no_warning = false; } - } // c - } // b + } + + } } // a } // ia2 } // it2 } // ia1 } // it1 + delete[] latvec; ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); if (!all_pass || !no_warning) { std::stringstream mess; From 442e95c43b19a3cfee142408cc39d13e7a29b31b Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Mon, 3 Mar 2025 14:40:53 +0800 Subject: [PATCH 03/14] add the change --- source/module_cell/check_atomic_stru.cpp | 94 ++++++++++++++++++------ 1 file changed, 73 insertions(+), 21 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index d4c679db23..cdf6116210 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -2,6 +2,10 @@ #include "module_base/timer.h" #include "module_base/element_covalent_radius.h" +void process_i(int i, double dx_lat, double dy_lat, double dz_lat) +{ + +} void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { // First we calculate all bond length in the structure, // and compare with the covalent_bond_length, @@ -48,7 +52,9 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) latvec[0]=ucell.a1.x; latvec[1]=ucell.a2.x; latvec[2]=ucell.a3.x; latvec[3]=ucell.a1.y; latvec[4]=ucell.a2.y; latvec[5]=ucell.a3.y; latvec[6]=ucell.a1.z; latvec[7]=ucell.a2.z; latvec[8]=ucell.a3.z; - std::vector A(27*3); + double* A = new double[27*3]; + std::vector cell(27); + std::vector label(ntype); for (int i=0;i<27;i++) { int a = (i / 9) % 3 - 1; // 计算a的值 @@ -57,6 +63,17 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) A[3*i]= a*latvec[0] + b *latvec[1] + c* latvec[2]; A[3*i+1]= a*latvec[3] + b *latvec[4] + c* latvec[5]; A[3*i+2]= a*latvec[6] + b *latvec[7] + c* latvec[8]; + std::ostringstream tmp_oss; + tmp_oss << " (cell:" << std::setw(2) << a + << " " << std::setw(2) << b << " " + << std::setw(2) << c << "), distance= "; + cell[i] = tmp_oss.str(); + } + for (int it=0;it Date: Mon, 3 Mar 2025 19:58:47 +0800 Subject: [PATCH 04/14] change --- source/driver_run.cpp | 78 +++--- source/module_cell/CMakeLists.txt | 1 + source/module_cell/check_atomic_stru.cpp | 237 +++++++--------- source/module_cell/check_atomic_stru.h | 4 + source/module_cell/check_atomic_stru1.cpp | 317 ++++++++++++++++++++++ 5 files changed, 459 insertions(+), 178 deletions(-) create mode 100644 source/module_cell/check_atomic_stru1.cpp diff --git a/source/driver_run.cpp b/source/driver_run.cpp index 9955503d90..e3401f820f 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -48,51 +48,51 @@ void Driver::driver_run() ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); - //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) - ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); +// //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) +// ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); - //! 3: initialize Esolver and fill json-structure - p_esolver->before_all_runners(ucell, PARAM.inp); +// //! 3: initialize Esolver and fill json-structure +// p_esolver->before_all_runners(ucell, PARAM.inp); - // this Json part should be moved to before_all_runners, mohan 2024-05-12 -#ifdef __RAPIDJSON - Json::gen_stru_wrapper(&ucell); -#endif +// // this Json part should be moved to before_all_runners, mohan 2024-05-12 +// #ifdef __RAPIDJSON +// Json::gen_stru_wrapper(&ucell); +// #endif - const std::string cal_type = PARAM.inp.calculation; +// const std::string cal_type = PARAM.inp.calculation; - //! 4: different types of calculations - if (cal_type == "md") - { - Run_MD::md_line(ucell, p_esolver, PARAM); - } - else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") - { - Relax_Driver rl_driver; - rl_driver.relax_driver(p_esolver, ucell); - } - else if (cal_type == "get_S") - { - p_esolver->runner(ucell, 0); - } - else - { - //! supported "other" functions: - //! get_pchg(LCAO), - //! test_memory(PW,LCAO), - //! test_neighbour(LCAO), - //! gen_bessel(PW), et al. - const int istep = 0; - p_esolver->others(ucell, istep); - } +// //! 4: different types of calculations +// if (cal_type == "md") +// { +// Run_MD::md_line(ucell, p_esolver, PARAM); +// } +// else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") +// { +// Relax_Driver rl_driver; +// rl_driver.relax_driver(p_esolver, ucell); +// } +// else if (cal_type == "get_S") +// { +// p_esolver->runner(ucell, 0); +// } +// else +// { +// //! supported "other" functions: +// //! get_pchg(LCAO), +// //! test_memory(PW,LCAO), +// //! test_neighbour(LCAO), +// //! gen_bessel(PW), et al. +// const int istep = 0; +// p_esolver->others(ucell, istep); +// } - //! 5: clean up esolver - p_esolver->after_all_runners(ucell); +// //! 5: clean up esolver +// p_esolver->after_all_runners(ucell); - ModuleESolver::clean_esolver(p_esolver); +// ModuleESolver::clean_esolver(p_esolver); - //! 6: output the json file - Json::create_Json(&ucell, PARAM); - ModuleBase::timer::tick("Driver", "driver_line"); +// //! 6: output the json file +// Json::create_Json(&ucell, PARAM); +// ModuleBase::timer::tick("Driver", "driver_line"); return; } diff --git a/source/module_cell/CMakeLists.txt b/source/module_cell/CMakeLists.txt index 3ce2763b5d..3e2e49e13e 100644 --- a/source/module_cell/CMakeLists.txt +++ b/source/module_cell/CMakeLists.txt @@ -23,6 +23,7 @@ add_library( parallel_kpoints.cpp cell_index.cpp check_atomic_stru.cpp + check_atomic_stru1.cpp update_cell.cpp bcast_cell.cpp read_stru.cpp diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index cdf6116210..f7f59c9806 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -1,12 +1,15 @@ #include "check_atomic_stru.h" -#include "module_base/timer.h" + #include "module_base/element_covalent_radius.h" +#include "module_base/timer.h" void process_i(int i, double dx_lat, double dy_lat, double dz_lat) { - } -void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { + + +void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) +{ // First we calculate all bond length in the structure, // and compare with the covalent_bond_length, // if there has bond length is shorter than covalent_bond_length * factor, @@ -18,13 +21,13 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) bool all_pass = true; bool no_warning = true; const double warning_coef = 0.6; - const double min_factor_coef=std::min(warning_coef,factor); + const double min_factor_coef = std::min(warning_coef, factor); std::stringstream errorlog; errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); std::vector symbol_covalent_radiuss(ntype); - for (int it = 0;it < ntype ; it++) + for (int it = 0; it < ntype; it++) { std::string symbol1 = ""; for (char ch: ucell.atoms[it].label) @@ -37,7 +40,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) { - symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); + symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); } else { @@ -47,160 +50,115 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) GlobalV::ofs_running << mess.str(); std::cout << mess.str(); } - } - double * latvec = new double [9]; - latvec[0]=ucell.a1.x; latvec[1]=ucell.a2.x; latvec[2]=ucell.a3.x; - latvec[3]=ucell.a1.y; latvec[4]=ucell.a2.y; latvec[5]=ucell.a3.y; - latvec[6]=ucell.a1.z; latvec[7]=ucell.a2.z; latvec[8]=ucell.a3.z; - double* A = new double[27*3]; + } + double* latvec = new double[9]; + latvec[0] = ucell.a1.x; + latvec[1] = ucell.a2.x; + latvec[2] = ucell.a3.x; + latvec[3] = ucell.a1.y; + latvec[4] = ucell.a2.y; + latvec[5] = ucell.a3.y; + latvec[6] = ucell.a1.z; + latvec[7] = ucell.a2.z; + latvec[8] = ucell.a3.z; + double* A = new double[27 * 3]; std::vector cell(27); std::vector label(ntype); - for (int i=0;i<27;i++) + for (int i = 0; i < 27; i++) { int a = (i / 9) % 3 - 1; // 计算a的值 int b = (i / 3) % 3 - 1; // 计算b的值 int c = i % 3 - 1; // 计算c的值 - A[3*i]= a*latvec[0] + b *latvec[1] + c* latvec[2]; - A[3*i+1]= a*latvec[3] + b *latvec[4] + c* latvec[5]; - A[3*i+2]= a*latvec[6] + b *latvec[7] + c* latvec[8]; + A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; + A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; + A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; std::ostringstream tmp_oss; - tmp_oss << " (cell:" << std::setw(2) << a - << " " << std::setw(2) << b << " " - << std::setw(2) << c << "), distance= "; + tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c + << "), distance= "; cell[i] = tmp_oss.str(); } - for (int it=0;it delta_lat(3); + // #pragma omp parallel for + for(int iat = 0; iat < ucell.nat; iat++) + { + const int it1 = ucell.iat2it[iat]; + const int ia1 = ucell.iat2ia[iat]; + const double symbol1_covalent_radius = symbol_covalent_radiuss[it1]; + double x1 = ucell.atoms[it1].taud[ia1].x; + double y1 = ucell.atoms[it1].taud[ia1].y; + double z1 = ucell.atoms[it1].taud[ia1].z; + for (int it2 = it1; it2 < ntype; it2++) + { + double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; + const double bohr_to_a= ModuleBase::BOHR_TO_A; + double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; + const double min_error = covalent_length * min_factor_coef / ucell.lat0; + const double min_error_2 = min_error * min_error; + const double factor_error = covalent_length * factor; + for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) + { + const bool is_same_atom = (it1 == it2) && (ia1 == ia2); + double delta_x = ucell.atoms[it2].taud[ia2].x - x1; + double delta_y = ucell.atoms[it2].taud[ia2].y - y1; + double delta_z = ucell.atoms[it2].taud[ia2].z - z1; + delta_lat[0] = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2]; + delta_lat[1] = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5]; + delta_lat[2] = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8]; + for (int i = 0; i < 27; i++) { - double delta_x = ucell.atoms[it2].taud[ia2].x - x1; - double delta_y = ucell.atoms[it2].taud[ia2].y - y1; - double delta_z = ucell.atoms[it2].taud[ia2].z - z1; - const double delta_x_lat = delta_x*latvec[0] + delta_y*latvec[1]+ delta_z*latvec[2]; - const double delta_y_lat = delta_x*latvec[3] + delta_y*latvec[4]+ delta_z*latvec[5]; - const double delta_z_lat = delta_x*latvec[6] + delta_y*latvec[7]+ delta_z*latvec[8]; - - double* current_ptr = &A[0]; - const bool is_same_atom = (it1 == it2) && (ia1 == ia2); - const bool warning_error = all_pass || no_warning; - if (is_same_atom) + if ((is_same_atom) && (i==13)) + continue; + const int offset = i*3; + const double part1 = delta_lat[0] + A[offset] ; + const double part2 = delta_lat[1] + A[offset+1]; + const double part3 = delta_lat[2] + A[offset+2]; + const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; + const bool flag = bond_length < min_error_2 ? true : false; + if (flag) { - - for (int i=0;i<27;i++) - { - if (i==13) - continue; - double part1= delta_x_lat+A[3*i]; - double part2= delta_y_lat+A[3*i+1]; - double part3= delta_z_lat+A[3*i+2]; - double bond_length = part1*part1 + part2*part2+ part3*part3; - - if (bond_length < min_error_2) + const double sqrt_bon = sqrt(bond_length) / lat0; + // #pragma omp critical { - const double sqrt_bon = sqrt(bond_length)/lat0; - errorlog << std::setw(3) << ia1 + 1 << "-th " - << label[it1] << ", " - << std::setw(3) << ia2 + 1 << "-th " - << label[it2] - << cell[i] - << std::setprecision(3) - << sqrt_bon << " Bohr (" - << sqrt_bon * ModuleBase::BOHR_TO_A - << " Angstrom)\n"; - if (all_pass || no_warning) - { - if (bond_length < pow(covalent_length * factor/lat0,2)) - { - all_pass = false; - } else { - no_warning = false; - } - } - + no_warning = false; + all_pass = bond_length < factor_error ? false :true; + errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) + << ia2 + 1 << "-th " << &label[it2] << &cell[i] + << std::setprecision(3) << sqrt_bon + << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; } - } - } - else - { - for (int i=0;i<27;i++) - { - const double part1= delta_x_lat + *(current_ptr++); - const double part2= delta_y_lat + *(current_ptr++); - const double part3= delta_z_lat + *(current_ptr++); - const double bond_length = part1*part1 + part2*part2+ part3*part3; - - if (bond_length < min_error_2) - { - const double sqrt_bon = sqrt(bond_length)/lat0; - errorlog << std::setw(3) << ia1 + 1 << "-th " - << label[it1] << ", " - << std::setw(3) << ia2 + 1 << "-th " - << label[it2] - << cell[i] - << std::setprecision(3) - << sqrt_bon << " Bohr (" - << sqrt_bon * ModuleBase::BOHR_TO_A - << " Angstrom)\n"; - if (warning_error) - { - if (bond_length < pow(covalent_length * factor/lat0,2)) - { - all_pass = false; - } else { - no_warning = false; - } - } - - } - } // a } - - } // ia2 - } // it2 - } // ia1 - } // it1 + + } // a + + } // ia2 + } // it2 + } delete[] latvec; ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - if (!all_pass || !no_warning) { + if (!all_pass || !no_warning) + { std::stringstream mess; - mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" - << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - << std::endl; + mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; - mess << "!!! Please check the nearest-neighbor list in log file." - << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" - << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - << std::endl; + mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_running << mess.str() << mess.str() << mess.str() - << errorlog.str(); + GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); std::cout << mess.str() << mess.str() << mess.str() << std::endl; std::ofstream outFile("example2.txt"); - if (outFile.is_open()) { + if (outFile.is_open()) + { // 将stringstream中的内容写入到文件 outFile << errorlog.str(); @@ -208,21 +166,22 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) // 关闭文件 outFile.close(); std::cout << "数据成功写入到文件。\n"; - } else { + } + else + { std::cerr << "无法打开文件进行写入。\n"; } - if (!all_pass) { + if (!all_pass) + { mess.clear(); mess.str(""); mess << "If this structure is what you want, you can set " "'min_dist_coef'" << std::endl; - mess << "as a smaller value (the current value is " << factor - << ") in INPUT file." << std::endl; + mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; GlobalV::ofs_running << mess.str(); std::cout << mess.str(); // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } } - } diff --git a/source/module_cell/check_atomic_stru.h b/source/module_cell/check_atomic_stru.h index 03e69918d4..145e0599eb 100644 --- a/source/module_cell/check_atomic_stru.h +++ b/source/module_cell/check_atomic_stru.h @@ -6,6 +6,10 @@ class Check_Atomic_Stru { public: static void check_atomic_stru(UnitCell& ucell, const double& factor); + + static void check_atomic_stru1(UnitCell& ucell, const double& factor); + + static void check_atomic_stru2(UnitCell& ucell, const double& factor); }; #endif \ No newline at end of file diff --git a/source/module_cell/check_atomic_stru1.cpp b/source/module_cell/check_atomic_stru1.cpp new file mode 100644 index 0000000000..8383220529 --- /dev/null +++ b/source/module_cell/check_atomic_stru1.cpp @@ -0,0 +1,317 @@ +#include "check_atomic_stru.h" +#include "module_base/timer.h" +#include "module_base/element_covalent_radius.h" +#include "module_base/timer.h" +#include "omp.h" +void Check_Atomic_Stru::check_atomic_stru2(UnitCell& ucell, const double& factor) +{ + // First we calculate all bond length in the structure, + // and compare with the covalent_bond_length, + // if there has bond length is shorter than covalent_bond_length * factor, + // we think this structure is unreasonable. + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); + assert(ucell.ntype > 0); + const int ntype = ucell.ntype; + const double lat0 = ucell.lat0; + bool all_pass = true; + bool no_warning = true; + const double warning_coef = 0.6; + const double min_factor_coef = std::min(warning_coef, factor); + + std::stringstream errorlog; + errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); + + std::vector symbol_covalent_radiuss(ntype); + for (int it = 0; it < ntype; it++) + { + std::string symbol1 = ""; + for (char ch: ucell.atoms[it].label) + { + if (std::isalpha(ch)) + { + symbol1.push_back(ch); + } + } + + if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) + { + symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); + } + else + { + std::stringstream mess; + mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; + mess << "set the covalent radius to be 0." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + } + } + double* latvec = new double[9]; + latvec[0] = ucell.a1.x; + latvec[1] = ucell.a2.x; + latvec[2] = ucell.a3.x; + latvec[3] = ucell.a1.y; + latvec[4] = ucell.a2.y; + latvec[5] = ucell.a3.y; + latvec[6] = ucell.a1.z; + latvec[7] = ucell.a2.z; + latvec[8] = ucell.a3.z; + double* A = new double[27 * 3]; + std::vector cell(27); + std::vector label(ntype); + for (int i = 0; i < 27; i++) + { + int a = (i / 9) % 3 - 1; // 计算a的值 + int b = (i / 3) % 3 - 1; // 计算b的值 + int c = i % 3 - 1; // 计算c的值 + A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; + A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; + A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; + std::ostringstream tmp_oss; + tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c + << "), distance= "; + cell[i] = tmp_oss.str(); + } + for (int it = 0; it < ntype; it++) + { + std::ostringstream tmp_oss; + tmp_oss << std::setw(3) << ucell.atoms[it].label; + label[it] = tmp_oss.str(); + } + #pragma omp parallel + { + #pragma omp single + for (int it1 = 0; it1 < ntype; it1++) + { + #pragma omp task firstprivate(it1) + for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) + { + const double symbol1_covalent_radius = symbol_covalent_radiuss[it1]; + double x1 = ucell.atoms[it1].taud[ia1].x; + double y1 = ucell.atoms[it1].taud[ia1].y; + double z1 = ucell.atoms[it1].taud[ia1].z; + + for (int it2 = it1; it2 < ntype; it2++) + { + double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; + const double bohr_to_a= ModuleBase::BOHR_TO_A; + double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; + const double min_error = covalent_length * min_factor_coef / ucell.lat0; + const double min_error_2 = min_error * min_error; + const double factor_error = covalent_length * factor / ucell.lat0; + for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) + { + double delta_x = ucell.atoms[it2].taud[ia2].x - x1; + double delta_y = ucell.atoms[it2].taud[ia2].y - y1; + double delta_z = ucell.atoms[it2].taud[ia2].z - z1; + const double delta_x_lat = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2]; + const double delta_y_lat = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5]; + const double delta_z_lat = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8]; + + double* current_ptr = &A[0]; + const bool is_same_atom = (it1 == it2) && (ia1 == ia2); + + for (int i = 0; i < 27; i++) + { + if (is_same_atom && i==13) + continue; + const double part1 = delta_x_lat + *(current_ptr++); + const double part2 = delta_y_lat + *(current_ptr++); + const double part3 = delta_z_lat + *(current_ptr++); + const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; + if (bond_length < min_error_2) + { + const double sqrt_bon = sqrt(bond_length) / lat0; + no_warning = false; + all_pass = bond_length < factor_error ? false :true; + // errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) + // << ia2 + 1 << "-th " << &label[it2] << &cell[i] + // << std::setprecision(3) << sqrt_bon + // << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; + } + } // a + + } // ia2 + } // it2 + } // ia1 + } // it1 + } + delete[] latvec; + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); + if (!all_pass || !no_warning) + { + std::stringstream mess; + mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; + mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + + GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); + std::cout << mess.str() << mess.str() << mess.str() << std::endl; + std::ofstream outFile("example2.txt"); + if (outFile.is_open()) + { + + // 将stringstream中的内容写入到文件 + outFile << errorlog.str(); + + // 关闭文件 + outFile.close(); + std::cout << "数据成功写入到文件。\n"; + } + else + { + std::cerr << "无法打开文件进行写入。\n"; + } + if (!all_pass) + { + mess.clear(); + mess.str(""); + mess << "If this structure is what you want, you can set " + "'min_dist_coef'" + << std::endl; + mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); + } + } +} + + +void matrix_multiply(const std::vector>& A, + const std::vector>& B, + std::vector>& C, + int N) { + // 使用OpenMP并行化外层循环，即对C中的每一行进行并行处理 + ModuleBase::timer::tick("Check_Atomic_Stru", "matrix_multiply"); + #pragma omp parallel for collapse(2) + for (int i = 0; i < N; ++i) { + for (int j = 0; j < N; ++j) { + double sum = 0.0; + for (int k = 0; k < N; ++k) { + sum += A[i][k] * B[k][j]; + } + C[i][j] = sum; + } + } + ModuleBase::timer::tick("Check_Atomic_Stru", "matrix_multiply"); +} + +void test() +{ + ModuleBase::timer::tick("Check_Atomic_Stru", "test"); + std::stringstream errorlog; + const double delta_x_lat=10; + const double delta_y_lat=10; + const double delta_z_lat=10; + const double min_error_2=1.0; + const double bohr_to_a=1.0; + const double lat0 =1.0; + double *A =new double [27*3]; + bool no_warning=0; + bool all_pass=0; + std::vector local_results(omp_get_num_threads()); + #pragma omp parallel + { + int tid=omp_get_thread_num(); + std::vector local_result; + #pragma omp parallel for + for (int j=0;j<1024*1024;j++) + for (int i = 0; i < 27; i++) + { + // if (is_same_atom && i==13) + // continue; + const double part1 = delta_x_lat + A[3*i]; + const double part2 = delta_y_lat + A[3*i+1]; + const double part3 = delta_z_lat + A[3*i+2]; + const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; + // if (bond_length < min_error_2) + // { + const double sqrt_bon = sqrt(bond_length) / lat0; + no_warning = false; + all_pass = bond_length < min_error_2 ? false :true; + // errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) + // << ia2 + 1 << "-th " << &label[it2] << &cell[i] + // local_results[tid] << std::setprecision(3) << sqrt_bon + // << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; + local_result.push_back( " Bohr ("); + + // } + } + } + ModuleBase::timer::tick("Check_Atomic_Stru", "test"); +} + + +void process_data() { + ModuleBase::timer::tick("Check_Atomic_Stru", "process_data"); + std::vector local_results; + const size_t total_iterations = 1024 * 1024 *1024; + const double lat0 = 1.0; + const double bohr_to_a = 0.529177; // 示例值 + const double min_error_2 = 1.0; // 示例值 + + std::vector results; + results.reserve(total_iterations); // 预分配空间以提高性能 + + #pragma omp parallel + { + if (omp_get_thread_num() == 0) { + local_results.clear(); // 初始化线程局部存储 + } + + #pragma omp barrier // 确保所有线程都已初始化它们的局部存储 + + #pragma omp for + for (size_t i = 0; i < total_iterations; ++i) { + const double delta_x_lat = sqrt(i + i); + const double delta_y_lat = sqrt(i * i); + const double delta_z_lat = sqrt(i * i + 10); + const double sqrt_bon = sqrt(delta_x_lat + delta_y_lat + delta_z_lat) / lat0; + + if (sqrt_bon < min_error_2) { + std::stringstream ss; + ss << std::setprecision(3) << sqrt_bon + << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; + local_results.push_back(ss.str()); + } + } + + // 合并线程局部存储到全局存储 + #pragma omp critical + { + results.insert(results.end(), local_results.begin(), local_results.end()); + } + } + + // 统一输出结果 + for (const auto& result : results) { + std::cerr << result; + } + ModuleBase::timer::tick("Check_Atomic_Stru", "process_data"); +} + +// 注意：这里假设 local_results 已经在某个地方定义了 +std::vector local_results; + +void Check_Atomic_Stru::check_atomic_stru1(UnitCell& ucell, const double& factor) +{ + const double lat0=1.0; + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru1"); + const int N = 1000; // 矩阵大小 + + // 初始化矩阵A和B + std::vector> A(N, std::vector(N)); + std::vector> B(N, std::vector(N)); + std::vector> C(N, std::vector(N)); + // 执行矩阵乘法 + matrix_multiply(A, B, C, N); + // test(); + process_data(); + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru1"); +} \ No newline at end of file From 471fa502649179ede55d5e7a27ceda42a17105c4 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Tue, 4 Mar 2025 17:20:04 +0800 Subject: [PATCH 05/14] change gtest --- source/driver_run.cpp | 78 ++-- source/module_cell/check_atomic_stru.cpp | 15 - source/module_cell/check_atomic_stru.h | 4 - .../module_cell/test/unitcell_test_readpp.cpp | 395 +++++++++--------- 4 files changed, 234 insertions(+), 258 deletions(-) diff --git a/source/driver_run.cpp b/source/driver_run.cpp index e3401f820f..9955503d90 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -48,51 +48,51 @@ void Driver::driver_run() ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running); Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef); -// //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) -// ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); + //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`) + ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell); -// //! 3: initialize Esolver and fill json-structure -// p_esolver->before_all_runners(ucell, PARAM.inp); + //! 3: initialize Esolver and fill json-structure + p_esolver->before_all_runners(ucell, PARAM.inp); -// // this Json part should be moved to before_all_runners, mohan 2024-05-12 -// #ifdef __RAPIDJSON -// Json::gen_stru_wrapper(&ucell); -// #endif + // this Json part should be moved to before_all_runners, mohan 2024-05-12 +#ifdef __RAPIDJSON + Json::gen_stru_wrapper(&ucell); +#endif -// const std::string cal_type = PARAM.inp.calculation; + const std::string cal_type = PARAM.inp.calculation; -// //! 4: different types of calculations -// if (cal_type == "md") -// { -// Run_MD::md_line(ucell, p_esolver, PARAM); -// } -// else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") -// { -// Relax_Driver rl_driver; -// rl_driver.relax_driver(p_esolver, ucell); -// } -// else if (cal_type == "get_S") -// { -// p_esolver->runner(ucell, 0); -// } -// else -// { -// //! supported "other" functions: -// //! get_pchg(LCAO), -// //! test_memory(PW,LCAO), -// //! test_neighbour(LCAO), -// //! gen_bessel(PW), et al. -// const int istep = 0; -// p_esolver->others(ucell, istep); -// } + //! 4: different types of calculations + if (cal_type == "md") + { + Run_MD::md_line(ucell, p_esolver, PARAM); + } + else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf") + { + Relax_Driver rl_driver; + rl_driver.relax_driver(p_esolver, ucell); + } + else if (cal_type == "get_S") + { + p_esolver->runner(ucell, 0); + } + else + { + //! supported "other" functions: + //! get_pchg(LCAO), + //! test_memory(PW,LCAO), + //! test_neighbour(LCAO), + //! gen_bessel(PW), et al. + const int istep = 0; + p_esolver->others(ucell, istep); + } -// //! 5: clean up esolver -// p_esolver->after_all_runners(ucell); + //! 5: clean up esolver + p_esolver->after_all_runners(ucell); -// ModuleESolver::clean_esolver(p_esolver); + ModuleESolver::clean_esolver(p_esolver); -// //! 6: output the json file -// Json::create_Json(&ucell, PARAM); -// ModuleBase::timer::tick("Driver", "driver_line"); + //! 6: output the json file + Json::create_Json(&ucell, PARAM); + ModuleBase::timer::tick("Driver", "driver_line"); return; } diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index f7f59c9806..3d1dbcafa9 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -156,21 +156,6 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); std::cout << mess.str() << mess.str() << mess.str() << std::endl; - std::ofstream outFile("example2.txt"); - if (outFile.is_open()) - { - - // 将stringstream中的内容写入到文件 - outFile << errorlog.str(); - - // 关闭文件 - outFile.close(); - std::cout << "数据成功写入到文件。\n"; - } - else - { - std::cerr << "无法打开文件进行写入。\n"; - } if (!all_pass) { mess.clear(); diff --git a/source/module_cell/check_atomic_stru.h b/source/module_cell/check_atomic_stru.h index 145e0599eb..03e69918d4 100644 --- a/source/module_cell/check_atomic_stru.h +++ b/source/module_cell/check_atomic_stru.h @@ -6,10 +6,6 @@ class Check_Atomic_Stru { public: static void check_atomic_stru(UnitCell& ucell, const double& factor); - - static void check_atomic_stru1(UnitCell& ucell, const double& factor); - - static void check_atomic_stru2(UnitCell& ucell, const double& factor); }; #endif \ No newline at end of file diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index 4e40f217e2..52a0adda74 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/source/module_cell/test/unitcell_test_readpp.cpp @@ -343,16 +343,14 @@ TEST_F(UcellDeathTest, CheckStructure) { // trial 1 testing::internal::CaptureStdout(); double factor = 0.2; + ucell->set_iat2itia(); EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, - testing::HasSubstr("WARNING: Some atoms are too close!!!")); + EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!")); // trial 2 testing::internal::CaptureStdout(); factor = 0.4; - EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor), - ::testing::ExitedWithCode(1), - ""); + EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor),::testing::ExitedWithCode(1),""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("The structure is unreasonable!")); // trial 3 @@ -361,207 +359,204 @@ TEST_F(UcellDeathTest, CheckStructure) { factor = 0.2; EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); output = testing::internal::GetCapturedStdout(); - EXPECT_THAT( - output, - testing::HasSubstr("Notice: symbol 'arbitrary' is not an element " - "symbol!!!! set the covalent radius to be 0.")); + EXPECT_THAT(output,testing::HasSubstr("Notice: symbol 'arbitrary' is not an element " + "symbol!!!! set the covalent radius to be 0.")); // trial 4 ucell->atoms[0].label = "Fe1"; testing::internal::CaptureStdout(); factor = 0.2; EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, - testing::HasSubstr("WARNING: Some atoms are too close!!!")); -} - -TEST_F(UcellDeathTest, ReadPseudoWarning1) { - PARAM.input.pseudo_dir = pp_dir; - PARAM.input.out_element_info = true; - ucell->pseudo_fn[1] = "H_sr_lda.upf"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::read_pseudo(ofs, *ucell), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, - testing::HasSubstr("All DFT functional must consistent.")); -} - -TEST_F(UcellDeathTest, ReadPseudoWarning2) { - PARAM.input.pseudo_dir = pp_dir; - PARAM.input.out_element_info = true; - ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; - testing::internal::CaptureStdout(); - EXPECT_NO_THROW(elecstate::read_pseudo(ofs, *ucell)); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT( - output, - testing::HasSubstr("Warning: the number of valence electrons in " - "pseudopotential > 3 for Al: [Ne] 3s2 3p1")); -} - -TEST_F(UcellTest, CalNelec) { - elecstate::read_cell_pseudopots(pp_dir, ofs, *ucell); - EXPECT_EQ(4, ucell->atoms[0].ncpp.zv); - EXPECT_EQ(1, ucell->atoms[1].ncpp.zv); - EXPECT_EQ(1, ucell->atoms[0].na); - EXPECT_EQ(2, ucell->atoms[1].na); - double nelec = 0; - elecstate::cal_nelec(ucell->atoms, ucell->ntype, nelec); - EXPECT_DOUBLE_EQ(6, nelec); -} - -TEST_F(UcellTest, CalNbands) -{ - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 6); -} - -TEST_F(UcellTest, CalNbandsFractionElec) -{ - PARAM.input.nelec = 9.5; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 6); -} - -TEST_F(UcellTest, CalNbandsSOC) -{ - PARAM.input.lspinorb = true; - PARAM.input.nbands = 0; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 20); -} - -TEST_F(UcellTest, CalNbandsSDFT) -{ - PARAM.input.esolver_type = "sdft"; - std::vector nelec_spin(2, 5.0); - EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); -} - -TEST_F(UcellTest, CalNbandsLCAO) -{ - PARAM.input.basis_type = "lcao"; - std::vector nelec_spin(2, 5.0); - EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); -} - -TEST_F(UcellTest, CalNbandsLCAOINPW) -{ - PARAM.input.basis_type = "lcao_in_pw"; - PARAM.sys.nlocal = PARAM.input.nbands - 1; - std::vector nelec_spin(2, 5.0); - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("NLOCAL < NBANDS")); -} - -TEST_F(UcellTest, CalNbandsWarning1) -{ - PARAM.input.nbands = PARAM.input.nelec / 2 - 1; - std::vector nelec_spin(2, 5.0); - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("Too few bands!")); -} - -TEST_F(UcellTest, CalNbandsWarning2) -{ - PARAM.input.nspin = 2; - PARAM.input.nupdown = 4.0; - std::vector nelec_spin(2); - nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; - nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("Too few spin up bands!")); -} - -TEST_F(UcellTest, CalNbandsWarning3) -{ - PARAM.input.nspin = 2; - PARAM.input.nupdown = -4.0; - std::vector nelec_spin(2); - nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; - nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("Too few spin down bands!")); -} - -TEST_F(UcellTest, CalNbandsSpin1) -{ - PARAM.input.nspin = 1; - PARAM.input.nbands = 0; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 15); -} - -TEST_F(UcellTest, CalNbandsSpin1LCAO) -{ - PARAM.input.nspin = 1; - PARAM.input.nbands = 0; - PARAM.input.basis_type = "lcao"; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 6); -} - -TEST_F(UcellTest, CalNbandsSpin4) -{ - PARAM.input.nspin = 4; - PARAM.input.nbands = 0; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 30); -} - -TEST_F(UcellTest, CalNbandsSpin4LCAO) -{ - PARAM.input.nspin = 4; - PARAM.input.nbands = 0; - PARAM.input.basis_type = "lcao"; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 6); -} - -TEST_F(UcellTest, CalNbandsSpin2) -{ - PARAM.input.nspin = 2; - PARAM.input.nbands = 0; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 16); -} - -TEST_F(UcellTest, CalNbandsSpin2LCAO) -{ - PARAM.input.nspin = 2; - PARAM.input.nbands = 0; - PARAM.input.basis_type = "lcao"; - std::vector nelec_spin(2, 5.0); - elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); - EXPECT_EQ(PARAM.input.nbands, 6); -} - -TEST_F(UcellTest, CalNbandsGaussWarning) -{ - PARAM.input.nbands = 5; - std::vector nelec_spin(2, 5.0); - PARAM.input.smearing_method = "gaussian"; - testing::internal::CaptureStdout(); - EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); - EXPECT_THAT(output, testing::HasSubstr("for smearing, num. of bands > num. of occupied bands")); -} + EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!")); +} + +// TEST_F(UcellDeathTest, ReadPseudoWarning1) { +// PARAM.input.pseudo_dir = pp_dir; +// PARAM.input.out_element_info = true; +// ucell->pseudo_fn[1] = "H_sr_lda.upf"; +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::read_pseudo(ofs, *ucell), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, +// testing::HasSubstr("All DFT functional must consistent.")); +// } + +// TEST_F(UcellDeathTest, ReadPseudoWarning2) { +// PARAM.input.pseudo_dir = pp_dir; +// PARAM.input.out_element_info = true; +// ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; +// testing::internal::CaptureStdout(); +// EXPECT_NO_THROW(elecstate::read_pseudo(ofs, *ucell)); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT( +// output, +// testing::HasSubstr("Warning: the number of valence electrons in " +// "pseudopotential > 3 for Al: [Ne] 3s2 3p1")); +// } + +// TEST_F(UcellTest, CalNelec) { +// elecstate::read_cell_pseudopots(pp_dir, ofs, *ucell); +// EXPECT_EQ(4, ucell->atoms[0].ncpp.zv); +// EXPECT_EQ(1, ucell->atoms[1].ncpp.zv); +// EXPECT_EQ(1, ucell->atoms[0].na); +// EXPECT_EQ(2, ucell->atoms[1].na); +// double nelec = 0; +// elecstate::cal_nelec(ucell->atoms, ucell->ntype, nelec); +// EXPECT_DOUBLE_EQ(6, nelec); +// } + +// TEST_F(UcellTest, CalNbands) +// { +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 6); +// } + +// TEST_F(UcellTest, CalNbandsFractionElec) +// { +// PARAM.input.nelec = 9.5; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 6); +// } + +// TEST_F(UcellTest, CalNbandsSOC) +// { +// PARAM.input.lspinorb = true; +// PARAM.input.nbands = 0; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 20); +// } + +// TEST_F(UcellTest, CalNbandsSDFT) +// { +// PARAM.input.esolver_type = "sdft"; +// std::vector nelec_spin(2, 5.0); +// EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); +// } + +// TEST_F(UcellTest, CalNbandsLCAO) +// { +// PARAM.input.basis_type = "lcao"; +// std::vector nelec_spin(2, 5.0); +// EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); +// } + +// TEST_F(UcellTest, CalNbandsLCAOINPW) +// { +// PARAM.input.basis_type = "lcao_in_pw"; +// PARAM.sys.nlocal = PARAM.input.nbands - 1; +// std::vector nelec_spin(2, 5.0); +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, testing::HasSubstr("NLOCAL < NBANDS")); +// } + +// TEST_F(UcellTest, CalNbandsWarning1) +// { +// PARAM.input.nbands = PARAM.input.nelec / 2 - 1; +// std::vector nelec_spin(2, 5.0); +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, testing::HasSubstr("Too few bands!")); +// } + +// TEST_F(UcellTest, CalNbandsWarning2) +// { +// PARAM.input.nspin = 2; +// PARAM.input.nupdown = 4.0; +// std::vector nelec_spin(2); +// nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; +// nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, testing::HasSubstr("Too few spin up bands!")); +// } + +// TEST_F(UcellTest, CalNbandsWarning3) +// { +// PARAM.input.nspin = 2; +// PARAM.input.nupdown = -4.0; +// std::vector nelec_spin(2); +// nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; +// nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, testing::HasSubstr("Too few spin down bands!")); +// } + +// TEST_F(UcellTest, CalNbandsSpin1) +// { +// PARAM.input.nspin = 1; +// PARAM.input.nbands = 0; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 15); +// } + +// TEST_F(UcellTest, CalNbandsSpin1LCAO) +// { +// PARAM.input.nspin = 1; +// PARAM.input.nbands = 0; +// PARAM.input.basis_type = "lcao"; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 6); +// } + +// TEST_F(UcellTest, CalNbandsSpin4) +// { +// PARAM.input.nspin = 4; +// PARAM.input.nbands = 0; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 30); +// } + +// TEST_F(UcellTest, CalNbandsSpin4LCAO) +// { +// PARAM.input.nspin = 4; +// PARAM.input.nbands = 0; +// PARAM.input.basis_type = "lcao"; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 6); +// } + +// TEST_F(UcellTest, CalNbandsSpin2) +// { +// PARAM.input.nspin = 2; +// PARAM.input.nbands = 0; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 16); +// } + +// TEST_F(UcellTest, CalNbandsSpin2LCAO) +// { +// PARAM.input.nspin = 2; +// PARAM.input.nbands = 0; +// PARAM.input.basis_type = "lcao"; +// std::vector nelec_spin(2, 5.0); +// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); +// EXPECT_EQ(PARAM.input.nbands, 6); +// } + +// TEST_F(UcellTest, CalNbandsGaussWarning) +// { +// PARAM.input.nbands = 5; +// std::vector nelec_spin(2, 5.0); +// PARAM.input.smearing_method = "gaussian"; +// testing::internal::CaptureStdout(); +// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); +// output = testing::internal::GetCapturedStdout(); +// EXPECT_THAT(output, testing::HasSubstr("for smearing, num. of bands > num. of occupied bands")); +// } #ifdef __MPI #include "mpi.h" From 3b2e1f235f5d02eb235d686f31f928ee3e59da74 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Wed, 5 Mar 2025 14:40:18 +0800 Subject: [PATCH 06/14] add change --- source/module_cell/CMakeLists.txt | 1 - source/module_cell/check_atomic_stru.cpp | 5 +- source/module_cell/check_atomic_stru1.cpp | 317 ---------------------- 3 files changed, 2 insertions(+), 321 deletions(-) delete mode 100644 source/module_cell/check_atomic_stru1.cpp diff --git a/source/module_cell/CMakeLists.txt b/source/module_cell/CMakeLists.txt index 3e2e49e13e..3ce2763b5d 100644 --- a/source/module_cell/CMakeLists.txt +++ b/source/module_cell/CMakeLists.txt @@ -23,7 +23,6 @@ add_library( parallel_kpoints.cpp cell_index.cpp check_atomic_stru.cpp - check_atomic_stru1.cpp update_cell.cpp bcast_cell.cpp read_stru.cpp diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 3d1dbcafa9..ab6a286735 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -127,16 +127,14 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) // #pragma omp critical { no_warning = false; - all_pass = bond_length < factor_error ? false :true; + all_pass = all_pass && ( bond_length < factor_error ? false :true); errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) << ia2 + 1 << "-th " << &label[it2] << &cell[i] << std::setprecision(3) << sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; } } - } // a - } // ia2 } // it2 } @@ -169,4 +167,5 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } } + ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } diff --git a/source/module_cell/check_atomic_stru1.cpp b/source/module_cell/check_atomic_stru1.cpp deleted file mode 100644 index 8383220529..0000000000 --- a/source/module_cell/check_atomic_stru1.cpp +++ /dev/null @@ -1,317 +0,0 @@ -#include "check_atomic_stru.h" -#include "module_base/timer.h" -#include "module_base/element_covalent_radius.h" -#include "module_base/timer.h" -#include "omp.h" -void Check_Atomic_Stru::check_atomic_stru2(UnitCell& ucell, const double& factor) -{ - // First we calculate all bond length in the structure, - // and compare with the covalent_bond_length, - // if there has bond length is shorter than covalent_bond_length * factor, - // we think this structure is unreasonable. - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - assert(ucell.ntype > 0); - const int ntype = ucell.ntype; - const double lat0 = ucell.lat0; - bool all_pass = true; - bool no_warning = true; - const double warning_coef = 0.6; - const double min_factor_coef = std::min(warning_coef, factor); - - std::stringstream errorlog; - errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); - - std::vector symbol_covalent_radiuss(ntype); - for (int it = 0; it < ntype; it++) - { - std::string symbol1 = ""; - for (char ch: ucell.atoms[it].label) - { - if (std::isalpha(ch)) - { - symbol1.push_back(ch); - } - } - - if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) - { - symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); - } - else - { - std::stringstream mess; - mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; - mess << "set the covalent radius to be 0." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - } - } - double* latvec = new double[9]; - latvec[0] = ucell.a1.x; - latvec[1] = ucell.a2.x; - latvec[2] = ucell.a3.x; - latvec[3] = ucell.a1.y; - latvec[4] = ucell.a2.y; - latvec[5] = ucell.a3.y; - latvec[6] = ucell.a1.z; - latvec[7] = ucell.a2.z; - latvec[8] = ucell.a3.z; - double* A = new double[27 * 3]; - std::vector cell(27); - std::vector label(ntype); - for (int i = 0; i < 27; i++) - { - int a = (i / 9) % 3 - 1; // 计算a的值 - int b = (i / 3) % 3 - 1; // 计算b的值 - int c = i % 3 - 1; // 计算c的值 - A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; - A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; - A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; - std::ostringstream tmp_oss; - tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c - << "), distance= "; - cell[i] = tmp_oss.str(); - } - for (int it = 0; it < ntype; it++) - { - std::ostringstream tmp_oss; - tmp_oss << std::setw(3) << ucell.atoms[it].label; - label[it] = tmp_oss.str(); - } - #pragma omp parallel - { - #pragma omp single - for (int it1 = 0; it1 < ntype; it1++) - { - #pragma omp task firstprivate(it1) - for (int ia1 = 0; ia1 < ucell.atoms[it1].na; ia1++) - { - const double symbol1_covalent_radius = symbol_covalent_radiuss[it1]; - double x1 = ucell.atoms[it1].taud[ia1].x; - double y1 = ucell.atoms[it1].taud[ia1].y; - double z1 = ucell.atoms[it1].taud[ia1].z; - - for (int it2 = it1; it2 < ntype; it2++) - { - double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; - const double bohr_to_a= ModuleBase::BOHR_TO_A; - double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; - const double min_error = covalent_length * min_factor_coef / ucell.lat0; - const double min_error_2 = min_error * min_error; - const double factor_error = covalent_length * factor / ucell.lat0; - for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) - { - double delta_x = ucell.atoms[it2].taud[ia2].x - x1; - double delta_y = ucell.atoms[it2].taud[ia2].y - y1; - double delta_z = ucell.atoms[it2].taud[ia2].z - z1; - const double delta_x_lat = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2]; - const double delta_y_lat = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5]; - const double delta_z_lat = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8]; - - double* current_ptr = &A[0]; - const bool is_same_atom = (it1 == it2) && (ia1 == ia2); - - for (int i = 0; i < 27; i++) - { - if (is_same_atom && i==13) - continue; - const double part1 = delta_x_lat + *(current_ptr++); - const double part2 = delta_y_lat + *(current_ptr++); - const double part3 = delta_z_lat + *(current_ptr++); - const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; - if (bond_length < min_error_2) - { - const double sqrt_bon = sqrt(bond_length) / lat0; - no_warning = false; - all_pass = bond_length < factor_error ? false :true; - // errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) - // << ia2 + 1 << "-th " << &label[it2] << &cell[i] - // << std::setprecision(3) << sqrt_bon - // << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; - } - } // a - - } // ia2 - } // it2 - } // ia1 - } // it1 - } - delete[] latvec; - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - if (!all_pass || !no_warning) - { - std::stringstream mess; - mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; - mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - - GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); - std::cout << mess.str() << mess.str() << mess.str() << std::endl; - std::ofstream outFile("example2.txt"); - if (outFile.is_open()) - { - - // 将stringstream中的内容写入到文件 - outFile << errorlog.str(); - - // 关闭文件 - outFile.close(); - std::cout << "数据成功写入到文件。\n"; - } - else - { - std::cerr << "无法打开文件进行写入。\n"; - } - if (!all_pass) - { - mess.clear(); - mess.str(""); - mess << "If this structure is what you want, you can set " - "'min_dist_coef'" - << std::endl; - mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); - } - } -} - - -void matrix_multiply(const std::vector>& A, - const std::vector>& B, - std::vector>& C, - int N) { - // 使用OpenMP并行化外层循环，即对C中的每一行进行并行处理 - ModuleBase::timer::tick("Check_Atomic_Stru", "matrix_multiply"); - #pragma omp parallel for collapse(2) - for (int i = 0; i < N; ++i) { - for (int j = 0; j < N; ++j) { - double sum = 0.0; - for (int k = 0; k < N; ++k) { - sum += A[i][k] * B[k][j]; - } - C[i][j] = sum; - } - } - ModuleBase::timer::tick("Check_Atomic_Stru", "matrix_multiply"); -} - -void test() -{ - ModuleBase::timer::tick("Check_Atomic_Stru", "test"); - std::stringstream errorlog; - const double delta_x_lat=10; - const double delta_y_lat=10; - const double delta_z_lat=10; - const double min_error_2=1.0; - const double bohr_to_a=1.0; - const double lat0 =1.0; - double *A =new double [27*3]; - bool no_warning=0; - bool all_pass=0; - std::vector local_results(omp_get_num_threads()); - #pragma omp parallel - { - int tid=omp_get_thread_num(); - std::vector local_result; - #pragma omp parallel for - for (int j=0;j<1024*1024;j++) - for (int i = 0; i < 27; i++) - { - // if (is_same_atom && i==13) - // continue; - const double part1 = delta_x_lat + A[3*i]; - const double part2 = delta_y_lat + A[3*i+1]; - const double part3 = delta_z_lat + A[3*i+2]; - const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; - // if (bond_length < min_error_2) - // { - const double sqrt_bon = sqrt(bond_length) / lat0; - no_warning = false; - all_pass = bond_length < min_error_2 ? false :true; - // errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) - // << ia2 + 1 << "-th " << &label[it2] << &cell[i] - // local_results[tid] << std::setprecision(3) << sqrt_bon - // << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; - local_result.push_back( " Bohr ("); - - // } - } - } - ModuleBase::timer::tick("Check_Atomic_Stru", "test"); -} - - -void process_data() { - ModuleBase::timer::tick("Check_Atomic_Stru", "process_data"); - std::vector local_results; - const size_t total_iterations = 1024 * 1024 *1024; - const double lat0 = 1.0; - const double bohr_to_a = 0.529177; // 示例值 - const double min_error_2 = 1.0; // 示例值 - - std::vector results; - results.reserve(total_iterations); // 预分配空间以提高性能 - - #pragma omp parallel - { - if (omp_get_thread_num() == 0) { - local_results.clear(); // 初始化线程局部存储 - } - - #pragma omp barrier // 确保所有线程都已初始化它们的局部存储 - - #pragma omp for - for (size_t i = 0; i < total_iterations; ++i) { - const double delta_x_lat = sqrt(i + i); - const double delta_y_lat = sqrt(i * i); - const double delta_z_lat = sqrt(i * i + 10); - const double sqrt_bon = sqrt(delta_x_lat + delta_y_lat + delta_z_lat) / lat0; - - if (sqrt_bon < min_error_2) { - std::stringstream ss; - ss << std::setprecision(3) << sqrt_bon - << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; - local_results.push_back(ss.str()); - } - } - - // 合并线程局部存储到全局存储 - #pragma omp critical - { - results.insert(results.end(), local_results.begin(), local_results.end()); - } - } - - // 统一输出结果 - for (const auto& result : results) { - std::cerr << result; - } - ModuleBase::timer::tick("Check_Atomic_Stru", "process_data"); -} - -// 注意：这里假设 local_results 已经在某个地方定义了 -std::vector local_results; - -void Check_Atomic_Stru::check_atomic_stru1(UnitCell& ucell, const double& factor) -{ - const double lat0=1.0; - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru1"); - const int N = 1000; // 矩阵大小 - - // 初始化矩阵A和B - std::vector> A(N, std::vector(N)); - std::vector> B(N, std::vector(N)); - std::vector> C(N, std::vector(N)); - // 执行矩阵乘法 - matrix_multiply(A, B, C, N); - // test(); - process_data(); - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru1"); -} \ No newline at end of file From a407f06c5c22389bf6f621a209e8c7f378e1a55c Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Wed, 5 Mar 2025 16:16:44 +0800 Subject: [PATCH 07/14] change the file --- source/module_cell/check_atomic_stru.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index ab6a286735..f598325ac2 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -21,7 +21,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) bool all_pass = true; bool no_warning = true; const double warning_coef = 0.6; - const double min_factor_coef = std::min(warning_coef, factor); + const double max_factor_coef = std::max(warning_coef, factor); std::stringstream errorlog; errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); @@ -66,9 +66,9 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) std::vector label(ntype); for (int i = 0; i < 27; i++) { - int a = (i / 9) % 3 - 1; // 计算a的值 - int b = (i / 3) % 3 - 1; // 计算b的值 - int c = i % 3 - 1; // 计算c的值 + int a = (i / 9) % 3 - 1; + int b = (i / 3) % 3 - 1; + int c = i % 3 - 1; A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; @@ -85,7 +85,8 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) } std::vector delta_lat(3); - // #pragma omp parallel for + const double bohr_to_a= ModuleBase::BOHR_TO_A; + #pragma omp parallel for for(int iat = 0; iat < ucell.nat; iat++) { const int it1 = ucell.iat2it[iat]; @@ -97,10 +98,9 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) for (int it2 = it1; it2 < ntype; it2++) { double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; - const double bohr_to_a= ModuleBase::BOHR_TO_A; double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; - const double min_error = covalent_length * min_factor_coef / ucell.lat0; - const double min_error_2 = min_error * min_error; + const double max_error = covalent_length * max_factor_coef / ucell.lat0; + const double max_error_2 = max_error * max_error; const double factor_error = covalent_length * factor; for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) { @@ -120,7 +120,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) const double part2 = delta_lat[1] + A[offset+1]; const double part3 = delta_lat[2] + A[offset+2]; const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; - const bool flag = bond_length < min_error_2 ? true : false; + const bool flag = bond_length < max_error_2 ? true : false; if (flag) { const double sqrt_bon = sqrt(bond_length) / lat0; From d268a3eef777fc5cfdb7c8cf96be00bf284669c9 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Wed, 5 Mar 2025 20:47:43 +0800 Subject: [PATCH 08/14] add format --- source/module_cell/check_atomic_stru.cpp | 278 ++++++++++++----------- 1 file changed, 141 insertions(+), 137 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index f598325ac2..15218d6b07 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -7,165 +7,169 @@ void process_i(int i, double dx_lat, double dy_lat, double dz_lat) { } - void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { // First we calculate all bond length in the structure, // and compare with the covalent_bond_length, // if there has bond length is shorter than covalent_bond_length * factor, // we think this structure is unreasonable. - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - assert(ucell.ntype > 0); - const int ntype = ucell.ntype; - const double lat0 = ucell.lat0; - bool all_pass = true; - bool no_warning = true; - const double warning_coef = 0.6; - const double max_factor_coef = std::max(warning_coef, factor); + if (GlobalV::MY_RANK == 0) + { + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); + assert(ucell.ntype > 0); + const int ntype = ucell.ntype; + const double lat0 = ucell.lat0; + bool all_pass = true; + bool no_warning = true; + const double warning_coef = 0.6; + const double max_factor_coef = std::max(warning_coef, factor); - std::stringstream errorlog; - errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); + std::stringstream errorlog; + errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); - std::vector symbol_covalent_radiuss(ntype); - for (int it = 0; it < ntype; it++) - { - std::string symbol1 = ""; - for (char ch: ucell.atoms[it].label) + std::vector symbol_covalent_radiuss(ntype); + for (int it = 0; it < ntype; it++) { - if (std::isalpha(ch)) + std::string symbol1 = ""; + for (char ch: ucell.atoms[it].label) { - symbol1.push_back(ch); + if (std::isalpha(ch)) + { + symbol1.push_back(ch); + } } - } - if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) + if (ModuleBase::CovalentRadius.find(symbol1) != ModuleBase::CovalentRadius.end()) + { + symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); + } + else + { + std::stringstream mess; + mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; + mess << "set the covalent radius to be 0." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + } + } + double* latvec = new double[9]; + latvec[0] = ucell.a1.x; + latvec[1] = ucell.a2.x; + latvec[2] = ucell.a3.x; + latvec[3] = ucell.a1.y; + latvec[4] = ucell.a2.y; + latvec[5] = ucell.a3.y; + latvec[6] = ucell.a1.z; + latvec[7] = ucell.a2.z; + latvec[8] = ucell.a3.z; + double* A = new double[27 * 3]; + std::vector cell(27); + std::vector label(ntype); + for (int i = 0; i < 27; i++) { - symbol_covalent_radiuss[it] = ModuleBase::CovalentRadius.at(symbol1); + int a = (i / 9) % 3 - 1; + int b = (i / 3) % 3 - 1; + int c = i % 3 - 1; + A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; + A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; + A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; + std::ostringstream tmp_oss; + tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c + << "), distance= "; + cell[i] = tmp_oss.str(); } - else + for (int it = 0; it < ntype; it++) { - std::stringstream mess; - mess << "Notice: symbol '" << symbol1 << "' is not an element symbol!!!! "; - mess << "set the covalent radius to be 0." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); + std::ostringstream tmp_oss; + tmp_oss << std::setw(3) << ucell.atoms[it].label; + label[it] = tmp_oss.str(); } - } - double* latvec = new double[9]; - latvec[0] = ucell.a1.x; - latvec[1] = ucell.a2.x; - latvec[2] = ucell.a3.x; - latvec[3] = ucell.a1.y; - latvec[4] = ucell.a2.y; - latvec[5] = ucell.a3.y; - latvec[6] = ucell.a1.z; - latvec[7] = ucell.a2.z; - latvec[8] = ucell.a3.z; - double* A = new double[27 * 3]; - std::vector cell(27); - std::vector label(ntype); - for (int i = 0; i < 27; i++) - { - int a = (i / 9) % 3 - 1; - int b = (i / 3) % 3 - 1; - int c = i % 3 - 1; - A[3 * i] = a * latvec[0] + b * latvec[1] + c * latvec[2]; - A[3 * i + 1] = a * latvec[3] + b * latvec[4] + c * latvec[5]; - A[3 * i + 2] = a * latvec[6] + b * latvec[7] + c * latvec[8]; - std::ostringstream tmp_oss; - tmp_oss << " (cell:" << std::setw(2) << a << " " << std::setw(2) << b << " " << std::setw(2) << c - << "), distance= "; - cell[i] = tmp_oss.str(); - } - for (int it = 0; it < ntype; it++) - { - std::ostringstream tmp_oss; - tmp_oss << std::setw(3) << ucell.atoms[it].label; - label[it] = tmp_oss.str(); - } - - std::vector delta_lat(3); - const double bohr_to_a= ModuleBase::BOHR_TO_A; - #pragma omp parallel for - for(int iat = 0; iat < ucell.nat; iat++) - { - const int it1 = ucell.iat2it[iat]; - const int ia1 = ucell.iat2ia[iat]; - const double symbol1_covalent_radius = symbol_covalent_radiuss[it1]; - double x1 = ucell.atoms[it1].taud[ia1].x; - double y1 = ucell.atoms[it1].taud[ia1].y; - double z1 = ucell.atoms[it1].taud[ia1].z; - for (int it2 = it1; it2 < ntype; it2++) + + const double bohr_to_a = ModuleBase::BOHR_TO_A; +#pragma omp parallel { - double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; - double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; - const double max_error = covalent_length * max_factor_coef / ucell.lat0; - const double max_error_2 = max_error * max_error; - const double factor_error = covalent_length * factor; - for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) + std::vector delta_lat(3); +#pragma omp for schedule(dynamic) + for (int iat = 0; iat < ucell.nat; iat++) { - const bool is_same_atom = (it1 == it2) && (ia1 == ia2); - double delta_x = ucell.atoms[it2].taud[ia2].x - x1; - double delta_y = ucell.atoms[it2].taud[ia2].y - y1; - double delta_z = ucell.atoms[it2].taud[ia2].z - z1; - delta_lat[0] = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2]; - delta_lat[1] = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5]; - delta_lat[2] = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8]; - for (int i = 0; i < 27; i++) + const int it1 = ucell.iat2it[iat]; + const int ia1 = ucell.iat2ia[iat]; + const double symbol1_covalent_radius = symbol_covalent_radiuss[it1]; + double x1 = ucell.atoms[it1].taud[ia1].x; + double y1 = ucell.atoms[it1].taud[ia1].y; + double z1 = ucell.atoms[it1].taud[ia1].z; + for (int it2 = it1; it2 < ntype; it2++) { - if ((is_same_atom) && (i==13)) - continue; - const int offset = i*3; - const double part1 = delta_lat[0] + A[offset] ; - const double part2 = delta_lat[1] + A[offset+1]; - const double part3 = delta_lat[2] + A[offset+2]; - const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; - const bool flag = bond_length < max_error_2 ? true : false; - if (flag) + double symbol2_covalent_radius = symbol_covalent_radiuss[it2]; + double covalent_length = (symbol1_covalent_radius + symbol2_covalent_radius) / bohr_to_a; + const double max_error = covalent_length * max_factor_coef / ucell.lat0; + const double max_error_2 = max_error * max_error; + const double factor_error = covalent_length * factor; + for (int ia2 = ia1; ia2 < ucell.atoms[it2].na; ia2++) { - const double sqrt_bon = sqrt(bond_length) / lat0; - // #pragma omp critical + const bool is_same_atom = (it1 == it2) && (ia1 == ia2); + double delta_x = ucell.atoms[it2].taud[ia2].x - x1; + double delta_y = ucell.atoms[it2].taud[ia2].y - y1; + double delta_z = ucell.atoms[it2].taud[ia2].z - z1; + delta_lat[0] = delta_x * latvec[0] + delta_y * latvec[1] + delta_z * latvec[2]; + delta_lat[1] = delta_x * latvec[3] + delta_y * latvec[4] + delta_z * latvec[5]; + delta_lat[2] = delta_x * latvec[6] + delta_y * latvec[7] + delta_z * latvec[8]; + for (int i = 0; i < 27; i++) { - no_warning = false; - all_pass = all_pass && ( bond_length < factor_error ? false :true); - errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) - << ia2 + 1 << "-th " << &label[it2] << &cell[i] - << std::setprecision(3) << sqrt_bon - << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; - } - } - } // a - } // ia2 - } // it2 - } - delete[] latvec; - ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - if (!all_pass || !no_warning) - { - std::stringstream mess; - mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; - mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - - GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); - std::cout << mess.str() << mess.str() << mess.str() << std::endl; - if (!all_pass) + if ((is_same_atom) && (i == 13)) + continue; + const int offset = i * 3; + const double part1 = delta_lat[0] + A[offset]; + const double part2 = delta_lat[1] + A[offset + 1]; + const double part3 = delta_lat[2] + A[offset + 2]; + const double bond_length = part1 * part1 + part2 * part2 + part3 * part3; + const bool flag = bond_length < max_error_2 ? true : false; + if (flag) + { + const double sqrt_bon = sqrt(bond_length) / lat0; + #pragma omp critical + { + no_warning = false; + all_pass = all_pass && (bond_length < factor_error ? false : true); + errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) + << ia2 + 1 << "-th " << &label[it2] << &cell[i] << std::setprecision(3) + << sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; + } + } + } // ia2 + } // it2 + } // iat + } + } + delete[] latvec; + delete[] A; + ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); + if (!all_pass || !no_warning) { - mess.clear(); - mess.str(""); - mess << "If this structure is what you want, you can set " - "'min_dist_coef'" - << std::endl; - mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - // ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); + std::stringstream mess; + mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; + mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + + GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); + std::cout << mess.str() << mess.str() << mess.str() << std::endl; + if (!all_pass) + { + mess.clear(); + mess.str(""); + mess << "If this structure is what you want, you can set " + "'min_dist_coef'" + << std::endl; + mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); + } } } - ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } From 2ad9f1e0cd9c0e8fa33ef0b8a8b30d8ff816d801 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Wed, 5 Mar 2025 20:54:16 +0800 Subject: [PATCH 09/14] modify back test --- .../module_cell/test/unitcell_test_readpp.cpp | 380 +++++++++--------- 1 file changed, 191 insertions(+), 189 deletions(-) diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index 52a0adda74..b71856e12a 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/source/module_cell/test/unitcell_test_readpp.cpp @@ -350,7 +350,9 @@ TEST_F(UcellDeathTest, CheckStructure) { // trial 2 testing::internal::CaptureStdout(); factor = 0.4; - EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor),::testing::ExitedWithCode(1),""); + EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor), + ::testing::ExitedWithCode(1), + ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("The structure is unreasonable!")); // trial 3 @@ -360,7 +362,7 @@ TEST_F(UcellDeathTest, CheckStructure) { EXPECT_NO_THROW(Check_Atomic_Stru::check_atomic_stru(*ucell, factor)); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output,testing::HasSubstr("Notice: symbol 'arbitrary' is not an element " - "symbol!!!! set the covalent radius to be 0.")); + "symbol!!!! set the covalent radius to be 0.")); // trial 4 ucell->atoms[0].label = "Fe1"; testing::internal::CaptureStdout(); @@ -370,193 +372,193 @@ TEST_F(UcellDeathTest, CheckStructure) { EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!")); } -// TEST_F(UcellDeathTest, ReadPseudoWarning1) { -// PARAM.input.pseudo_dir = pp_dir; -// PARAM.input.out_element_info = true; -// ucell->pseudo_fn[1] = "H_sr_lda.upf"; -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::read_pseudo(ofs, *ucell), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, -// testing::HasSubstr("All DFT functional must consistent.")); -// } - -// TEST_F(UcellDeathTest, ReadPseudoWarning2) { -// PARAM.input.pseudo_dir = pp_dir; -// PARAM.input.out_element_info = true; -// ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; -// testing::internal::CaptureStdout(); -// EXPECT_NO_THROW(elecstate::read_pseudo(ofs, *ucell)); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT( -// output, -// testing::HasSubstr("Warning: the number of valence electrons in " -// "pseudopotential > 3 for Al: [Ne] 3s2 3p1")); -// } - -// TEST_F(UcellTest, CalNelec) { -// elecstate::read_cell_pseudopots(pp_dir, ofs, *ucell); -// EXPECT_EQ(4, ucell->atoms[0].ncpp.zv); -// EXPECT_EQ(1, ucell->atoms[1].ncpp.zv); -// EXPECT_EQ(1, ucell->atoms[0].na); -// EXPECT_EQ(2, ucell->atoms[1].na); -// double nelec = 0; -// elecstate::cal_nelec(ucell->atoms, ucell->ntype, nelec); -// EXPECT_DOUBLE_EQ(6, nelec); -// } - -// TEST_F(UcellTest, CalNbands) -// { -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 6); -// } - -// TEST_F(UcellTest, CalNbandsFractionElec) -// { -// PARAM.input.nelec = 9.5; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 6); -// } - -// TEST_F(UcellTest, CalNbandsSOC) -// { -// PARAM.input.lspinorb = true; -// PARAM.input.nbands = 0; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 20); -// } - -// TEST_F(UcellTest, CalNbandsSDFT) -// { -// PARAM.input.esolver_type = "sdft"; -// std::vector nelec_spin(2, 5.0); -// EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); -// } - -// TEST_F(UcellTest, CalNbandsLCAO) -// { -// PARAM.input.basis_type = "lcao"; -// std::vector nelec_spin(2, 5.0); -// EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); -// } - -// TEST_F(UcellTest, CalNbandsLCAOINPW) -// { -// PARAM.input.basis_type = "lcao_in_pw"; -// PARAM.sys.nlocal = PARAM.input.nbands - 1; -// std::vector nelec_spin(2, 5.0); -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, testing::HasSubstr("NLOCAL < NBANDS")); -// } - -// TEST_F(UcellTest, CalNbandsWarning1) -// { -// PARAM.input.nbands = PARAM.input.nelec / 2 - 1; -// std::vector nelec_spin(2, 5.0); -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, testing::HasSubstr("Too few bands!")); -// } - -// TEST_F(UcellTest, CalNbandsWarning2) -// { -// PARAM.input.nspin = 2; -// PARAM.input.nupdown = 4.0; -// std::vector nelec_spin(2); -// nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; -// nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, testing::HasSubstr("Too few spin up bands!")); -// } - -// TEST_F(UcellTest, CalNbandsWarning3) -// { -// PARAM.input.nspin = 2; -// PARAM.input.nupdown = -4.0; -// std::vector nelec_spin(2); -// nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; -// nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, testing::HasSubstr("Too few spin down bands!")); -// } - -// TEST_F(UcellTest, CalNbandsSpin1) -// { -// PARAM.input.nspin = 1; -// PARAM.input.nbands = 0; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 15); -// } - -// TEST_F(UcellTest, CalNbandsSpin1LCAO) -// { -// PARAM.input.nspin = 1; -// PARAM.input.nbands = 0; -// PARAM.input.basis_type = "lcao"; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 6); -// } - -// TEST_F(UcellTest, CalNbandsSpin4) -// { -// PARAM.input.nspin = 4; -// PARAM.input.nbands = 0; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 30); -// } - -// TEST_F(UcellTest, CalNbandsSpin4LCAO) -// { -// PARAM.input.nspin = 4; -// PARAM.input.nbands = 0; -// PARAM.input.basis_type = "lcao"; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 6); -// } - -// TEST_F(UcellTest, CalNbandsSpin2) -// { -// PARAM.input.nspin = 2; -// PARAM.input.nbands = 0; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 16); -// } - -// TEST_F(UcellTest, CalNbandsSpin2LCAO) -// { -// PARAM.input.nspin = 2; -// PARAM.input.nbands = 0; -// PARAM.input.basis_type = "lcao"; -// std::vector nelec_spin(2, 5.0); -// elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); -// EXPECT_EQ(PARAM.input.nbands, 6); -// } - -// TEST_F(UcellTest, CalNbandsGaussWarning) -// { -// PARAM.input.nbands = 5; -// std::vector nelec_spin(2, 5.0); -// PARAM.input.smearing_method = "gaussian"; -// testing::internal::CaptureStdout(); -// EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); -// output = testing::internal::GetCapturedStdout(); -// EXPECT_THAT(output, testing::HasSubstr("for smearing, num. of bands > num. of occupied bands")); -// } +TEST_F(UcellDeathTest, ReadPseudoWarning1) { + PARAM.input.pseudo_dir = pp_dir; + PARAM.input.out_element_info = true; + ucell->pseudo_fn[1] = "H_sr_lda.upf"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::read_pseudo(ofs, *ucell), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, + testing::HasSubstr("All DFT functional must consistent.")); +} + +TEST_F(UcellDeathTest, ReadPseudoWarning2) { + PARAM.input.pseudo_dir = pp_dir; + PARAM.input.out_element_info = true; + ucell->pseudo_fn[0] = "Al_ONCV_PBE-1.0.upf"; + testing::internal::CaptureStdout(); + EXPECT_NO_THROW(elecstate::read_pseudo(ofs, *ucell)); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT( + output, + testing::HasSubstr("Warning: the number of valence electrons in " + "pseudopotential > 3 for Al: [Ne] 3s2 3p1")); +} + +TEST_F(UcellTest, CalNelec) { + elecstate::read_cell_pseudopots(pp_dir, ofs, *ucell); + EXPECT_EQ(4, ucell->atoms[0].ncpp.zv); + EXPECT_EQ(1, ucell->atoms[1].ncpp.zv); + EXPECT_EQ(1, ucell->atoms[0].na); + EXPECT_EQ(2, ucell->atoms[1].na); + double nelec = 0; + elecstate::cal_nelec(ucell->atoms, ucell->ntype, nelec); + EXPECT_DOUBLE_EQ(6, nelec); +} + +TEST_F(UcellTest, CalNbands) +{ + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 6); +} + +TEST_F(UcellTest, CalNbandsFractionElec) +{ + PARAM.input.nelec = 9.5; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 6); +} + +TEST_F(UcellTest, CalNbandsSOC) +{ + PARAM.input.lspinorb = true; + PARAM.input.nbands = 0; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 20); +} + +TEST_F(UcellTest, CalNbandsSDFT) +{ + PARAM.input.esolver_type = "sdft"; + std::vector nelec_spin(2, 5.0); + EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); +} + +TEST_F(UcellTest, CalNbandsLCAO) +{ + PARAM.input.basis_type = "lcao"; + std::vector nelec_spin(2, 5.0); + EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands)); +} + +TEST_F(UcellTest, CalNbandsLCAOINPW) +{ + PARAM.input.basis_type = "lcao_in_pw"; + PARAM.sys.nlocal = PARAM.input.nbands - 1; + std::vector nelec_spin(2, 5.0); + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("NLOCAL < NBANDS")); +} + +TEST_F(UcellTest, CalNbandsWarning1) +{ + PARAM.input.nbands = PARAM.input.nelec / 2 - 1; + std::vector nelec_spin(2, 5.0); + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("Too few bands!")); +} + +TEST_F(UcellTest, CalNbandsWarning2) +{ + PARAM.input.nspin = 2; + PARAM.input.nupdown = 4.0; + std::vector nelec_spin(2); + nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; + nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("Too few spin up bands!")); +} + +TEST_F(UcellTest, CalNbandsWarning3) +{ + PARAM.input.nspin = 2; + PARAM.input.nupdown = -4.0; + std::vector nelec_spin(2); + nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0; + nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0; + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("Too few spin down bands!")); +} + +TEST_F(UcellTest, CalNbandsSpin1) +{ + PARAM.input.nspin = 1; + PARAM.input.nbands = 0; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 15); +} + +TEST_F(UcellTest, CalNbandsSpin1LCAO) +{ + PARAM.input.nspin = 1; + PARAM.input.nbands = 0; + PARAM.input.basis_type = "lcao"; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 6); +} + +TEST_F(UcellTest, CalNbandsSpin4) +{ + PARAM.input.nspin = 4; + PARAM.input.nbands = 0; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 30); +} + +TEST_F(UcellTest, CalNbandsSpin4LCAO) +{ + PARAM.input.nspin = 4; + PARAM.input.nbands = 0; + PARAM.input.basis_type = "lcao"; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 6); +} + +TEST_F(UcellTest, CalNbandsSpin2) +{ + PARAM.input.nspin = 2; + PARAM.input.nbands = 0; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 16); +} + +TEST_F(UcellTest, CalNbandsSpin2LCAO) +{ + PARAM.input.nspin = 2; + PARAM.input.nbands = 0; + PARAM.input.basis_type = "lcao"; + std::vector nelec_spin(2, 5.0); + elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands); + EXPECT_EQ(PARAM.input.nbands, 6); +} + +TEST_F(UcellTest, CalNbandsGaussWarning) +{ + PARAM.input.nbands = 5; + std::vector nelec_spin(2, 5.0); + PARAM.input.smearing_method = "gaussian"; + testing::internal::CaptureStdout(); + EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), ""); + output = testing::internal::GetCapturedStdout(); + EXPECT_THAT(output, testing::HasSubstr("for smearing, num. of bands > num. of occupied bands")); +} #ifdef __MPI #include "mpi.h" From f5adb781f2f5b92d81bf99617434b564727ff768 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Thu, 6 Mar 2025 10:15:49 +0800 Subject: [PATCH 10/14] add the change in the check_atom_stru --- source/module_cell/check_atomic_stru.cpp | 78 ++++++++++++------------ 1 file changed, 38 insertions(+), 40 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 15218d6b07..50fa3514c7 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -3,30 +3,28 @@ #include "module_base/element_covalent_radius.h" #include "module_base/timer.h" -void process_i(int i, double dx_lat, double dy_lat, double dz_lat) -{ -} - void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { // First we calculate all bond length in the structure, // and compare with the covalent_bond_length, // if there has bond length is shorter than covalent_bond_length * factor, // we think this structure is unreasonable. + assert(ucell.ntype > 0); + bool all_pass = true; + bool no_warning = true; + std::stringstream errorlog; + errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); + if (GlobalV::MY_RANK == 0) { ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - assert(ucell.ntype > 0); + const int ntype = ucell.ntype; const double lat0 = ucell.lat0; - bool all_pass = true; - bool no_warning = true; + const double warning_coef = 0.6; const double max_factor_coef = std::max(warning_coef, factor); - std::stringstream errorlog; - errorlog.setf(std::ios_base::fixed, std::ios_base::floatfield); - std::vector symbol_covalent_radiuss(ntype); for (int it = 0; it < ntype; it++) { @@ -86,10 +84,10 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) } const double bohr_to_a = ModuleBase::BOHR_TO_A; -#pragma omp parallel +// #pragma omp parallel { std::vector delta_lat(3); -#pragma omp for schedule(dynamic) +// #pragma omp for schedule(dynamic) for (int iat = 0; iat < ucell.nat; iat++) { const int it1 = ucell.iat2it[iat]; @@ -126,11 +124,11 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) const bool flag = bond_length < max_error_2 ? true : false; if (flag) { - const double sqrt_bon = sqrt(bond_length) / lat0; - #pragma omp critical + const double sqrt_bon = sqrt(bond_length) * lat0; + // #pragma omp critical { no_warning = false; - all_pass = all_pass && (bond_length < factor_error ? false : true); + all_pass = all_pass && (sqrt_bon < factor_error ? false : true); errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) << ia2 + 1 << "-th " << &label[it2] << &cell[i] << std::setprecision(3) << sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; @@ -144,32 +142,32 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) delete[] latvec; delete[] A; ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); - if (!all_pass || !no_warning) - { - std::stringstream mess; - mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; - mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; - mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + } + if (!all_pass || !no_warning) + { + std::stringstream mess; + mess << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "!!! WARNING: Some atoms are too close!!!" << std::endl; + mess << "!!! Please check the nearest-neighbor list in log file." << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + mess << "%%%%%% WARNING WARNING WARNING WARNING WARNING %%%%%%" << std::endl; + mess << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); - std::cout << mess.str() << mess.str() << mess.str() << std::endl; - if (!all_pass) - { - mess.clear(); - mess.str(""); - mess << "If this structure is what you want, you can set " - "'min_dist_coef'" - << std::endl; - mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; - GlobalV::ofs_running << mess.str(); - std::cout << mess.str(); - ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); - } + GlobalV::ofs_running << mess.str() << mess.str() << mess.str() << errorlog.str(); + std::cout << mess.str() << mess.str() << mess.str() << std::endl; + if (!all_pass) + { + mess.clear(); + mess.str(""); + mess << "If this structure is what you want, you can set " + "'min_dist_coef'" + << std::endl; + mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; + GlobalV::ofs_running << mess.str(); + std::cout << mess.str(); + ModuleBase::WARNING_QUIT("Input", "The structure is unreasonable!"); } } } From 0dde41d601dc2a1d74d8060f374213d2e98aaaa4 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Thu, 6 Mar 2025 10:32:29 +0800 Subject: [PATCH 11/14] move openmp --- source/module_cell/check_atomic_stru.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 50fa3514c7..690190347a 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -21,7 +21,6 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) const int ntype = ucell.ntype; const double lat0 = ucell.lat0; - const double warning_coef = 0.6; const double max_factor_coef = std::max(warning_coef, factor); @@ -84,10 +83,10 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) } const double bohr_to_a = ModuleBase::BOHR_TO_A; -// #pragma omp parallel +#pragma omp parallel { std::vector delta_lat(3); -// #pragma omp for schedule(dynamic) +#pragma omp for schedule(dynamic) for (int iat = 0; iat < ucell.nat; iat++) { const int it1 = ucell.iat2it[iat]; @@ -125,7 +124,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) if (flag) { const double sqrt_bon = sqrt(bond_length) * lat0; - // #pragma omp critical + #pragma omp critical { no_warning = false; all_pass = all_pass && (sqrt_bon < factor_error ? false : true); From 110110478318c7917a7e7b0a7221bea79a952b8c Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Thu, 6 Mar 2025 17:04:46 +0800 Subject: [PATCH 12/14] add globlv for the question --- source/module_cell/check_atomic_stru.cpp | 8 +++----- source/module_cell/test/unitcell_test_readpp.cpp | 15 +++++++++++++-- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 690190347a..5080a0d31b 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -128,8 +128,8 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { no_warning = false; all_pass = all_pass && (sqrt_bon < factor_error ? false : true); - errorlog << std::setw(3) << ia1 + 1 << "-th " << &label[it1] << ", " << std::setw(3) - << ia2 + 1 << "-th " << &label[it2] << &cell[i] << std::setprecision(3) + errorlog << std::setw(3) << ia1 + 1 << "-th " << label[it1] << ", " << std::setw(3) + << ia2 + 1 << "-th " << label[it2] << cell[i] << std::setprecision(3) << sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; } } @@ -160,9 +160,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) { mess.clear(); mess.str(""); - mess << "If this structure is what you want, you can set " - "'min_dist_coef'" - << std::endl; + mess << "If this structure is what you want, you can set 'min_dist_coef'\n"; mess << "as a smaller value (the current value is " << factor << ") in INPUT file." << std::endl; GlobalV::ofs_running << mess.str(); std::cout << mess.str(); diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index b71856e12a..0f25123799 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/source/module_cell/test/unitcell_test_readpp.cpp @@ -11,6 +11,7 @@ #include "module_elecstate/read_pseudo.h" #include #include +#include "string.h" #ifdef __MPI #include "mpi.h" #endif @@ -341,6 +342,7 @@ TEST_F(UcellDeathTest, CheckStructure) { EXPECT_FALSE(ucell->atoms[0].ncpp.has_so); EXPECT_FALSE(ucell->atoms[1].ncpp.has_so); // trial 1 + testing::internal::CaptureStdout(); double factor = 0.2; ucell->set_iat2itia(); @@ -348,13 +350,22 @@ TEST_F(UcellDeathTest, CheckStructure) { output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output,testing::HasSubstr("WARNING: Some atoms are too close!!!")); // trial 2 - testing::internal::CaptureStdout(); + GlobalV::ofs_running.open("CheckStructure2.txt"); + ::testing::FLAGS_gtest_death_test_style = "threadsafe"; factor = 0.4; EXPECT_EXIT(Check_Atomic_Stru::check_atomic_stru(*ucell, factor), ::testing::ExitedWithCode(1), ""); - output = testing::internal::GetCapturedStdout(); + std::ifstream ifs("CheckStructure2.txt"); + if (ifs.is_open()) + { + std::string line; + while (std::getline(ifs, line)) { + output+=line; + } + } EXPECT_THAT(output, testing::HasSubstr("The structure is unreasonable!")); + GlobalV::ofs_running.open("running.log"); // trial 3 ucell->atoms[0].label = "arbitrary"; testing::internal::CaptureStdout(); From 4e9acb6ae22be5d56303679b10c122157c164c3c Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Fri, 7 Mar 2025 11:59:11 +0800 Subject: [PATCH 13/14] delete vector --- source/module_cell/check_atomic_stru.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index 5080a0d31b..c68da31b69 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -49,7 +49,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) std::cout << mess.str(); } } - double* latvec = new double[9]; + std::vector latvec (9); latvec[0] = ucell.a1.x; latvec[1] = ucell.a2.x; latvec[2] = ucell.a3.x; @@ -59,7 +59,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) latvec[6] = ucell.a1.z; latvec[7] = ucell.a2.z; latvec[8] = ucell.a3.z; - double* A = new double[27 * 3]; + std::vector A(27); std::vector cell(27); std::vector label(ntype); for (int i = 0; i < 27; i++) @@ -133,13 +133,11 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) << sqrt_bon << " Bohr (" << sqrt_bon * bohr_to_a << " Angstrom)\n"; } } - } // ia2 - } // it2 - } // iat - } + } + } // ia2 + } // it2 + } // iat } - delete[] latvec; - delete[] A; ModuleBase::timer::tick("Check_Atomic_Stru", "Check_Atomic_Stru"); } if (!all_pass || !no_warning) From 9a919137527028e5db3e518287e07804f5632d73 Mon Sep 17 00:00:00 2001 From: ubuntu <3158793232@qq.com> Date: Fri, 7 Mar 2025 12:56:38 +0800 Subject: [PATCH 14/14] update bug --- source/module_cell/check_atomic_stru.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/module_cell/check_atomic_stru.cpp b/source/module_cell/check_atomic_stru.cpp index c68da31b69..acdc30076f 100644 --- a/source/module_cell/check_atomic_stru.cpp +++ b/source/module_cell/check_atomic_stru.cpp @@ -59,7 +59,7 @@ void Check_Atomic_Stru::check_atomic_stru(UnitCell& ucell, const double& factor) latvec[6] = ucell.a1.z; latvec[7] = ucell.a2.z; latvec[8] = ucell.a3.z; - std::vector A(27); + std::vector A(27*3); std::vector cell(27); std::vector label(ntype); for (int i = 0; i < 27; i++)